/*==Evaluating Large-Scale Government Investments in Fertilizer Adoption: The Ethiopian Experience======*/
/*==============================Thomas Assefa, Ellen McCullough, Guush Berhane  ========================*/
	
	version 18.0
	
	*Preliminaries 
	cap log close
	estimates drop _all
	log using "Tables/Results_log", replace
	
	* Defining outcomes 
	loc outcome_all_datasets NPS_Dummy DAP_Dummy Fert_Dummy /// // Fertilizer adoption - binary outcome */  
		 w_App_rate w_DAP_rate w_Urea_rate w_NPS_rate w_N_rate w_P_rate  ///  //Application rates
		 CL_w_Yield_kgha_1 CL_w_Yield_kgha_2 CL_w_Yield_kgha_3  /// //Yield - crop level 
	
	loc outcome_area	CL_w_area_1 CL_w_area_2 CL_w_area_3 /// // Area planted - crop level 
						w_area_planted area_tot_TWM //Area planted total 
		
	loc outcome_FTF_only /// 
		CL_V_w_harvest_kg_1  CL_V_w_harvest_kg_2 CL_V_w_harvest_kg_3 ///  //Value of prodcution - crop level 
		TOT_V_harvest_kg  /// //Value of production total 
		DAP_Price_mr /// //DAP price 
		pcexp pcexp_aeu pcexpD_R paexpD_R //Consumption - HDDS and consumption expenditures  

	* Defining covariates  
	loc covars_all land_area_sqm_10thou Female Litrate
	loc covars_area  Female Litrate 	
	
	*Executing Estimation Procedures
	foreach data in FTF /*AgSS*/ LSMS {
		use `data'.dta, clear 
		keep if inlist(real(substr(string(A05, "%9.0f"), 1, 1)), 1, 3, 4, 7) //Restricting AgSS and LSMS datasets to the regiosn of the FTF data 
	
		* 
		forvalues c = 1/3 {
		cap drop Y`c'm vh`c'm
		g Y`c'm = CL_w_Yield_kgha_`c'==. 
		bysort id: egen Y`c'm_ever = max(Y`c'm)
		replace CL_w_Yield_kgha_`c' = . if Y`c'm_ever==1
		if "`data'" == "FTF" {
            g vh`c'm = CL_V_w_harvest_kg_`c' == .
            bysort id: egen vh`c'm_ever = max(vh`c'm)
            replace CL_V_w_harvest_kg_`c' = . if vh`c'm_ever == 1
        }
		}


	* For FTF dataset, include outcome_FTF_only
		if "`data'" == "FTF" {
			loc outcome `outcome_all_datasets'  `outcome_FTF_only'
		}
		else {
			loc outcome `outcome_all_datasets'
		}

	 * Loop through variables
		foreach v in `outcome' {
			*Define covariates to be used dynamically 
			if strpos("`outcome_area'", "`v'") {
				loc covariates Female Litrate 		
				}
			else {
				loc covariates `covars_all'
			}
				
	
			* Run csdid with for 90 and 95 % CI separately 
			csdid `v'  `covariates', ivar(id) time(year) gvar(first_treat) method(dripw) notyet wboot  agg(group)  rseed(19810310) cluster(A05) level(95)			
			estimates store `v'_`data'
			*Storing control group sds by cohort
			su `v' if (first_treat==0| first_treat==3) & year==1
			loc `v'_max_`data'_c1 = r(sd)
			display "`v'_max_`data'_c1: ``v'_max_`data'_c1'"
			loc `v'_min_`data'_c1 = ``v'_max_`data'_c1'*-1
			
			su `v' if first_treat==0 & year==2
			loc `v'_max_`data'_c2 = r(sd)
			display "`v'_max_`data'_c2: ``v'_max_`data'_c2'"
			loc `v'_min_`data'_c2 = ``v'_max_`data'_c2'*-1			
			
			* Store placeholder values for AgSS max/min SDs (cohort 1 and 2)
			loc `v'_max_AgSS_c1 = .
			display "`v'_max_AgSS_c1: ``v'_max_AgSS_c1'"
			loc `v'_min_AgSS_c1 = .

			loc `v'_max_AgSS_c2 = .
			display "`v'_max_AgSS_c2: ``v'_max_AgSS_c2'"
			loc `v'_min_AgSS_c2 = .

			*Storing 95% CIs 				
			matrix `v'_`data'_95 = e(cband)

			
			csdid `v'  `covariates', ivar(id) time(year) gvar(first_treat) method(dripw) notyet wboot  agg(group)  rseed(19810310) cluster(A05) level(90)			 
			*Storing 90% CIs 				
			matrix `v'_`data'_90 = e(cband)
					
			*Combining the 90 and 95 % CIs. 
			matrix `v'_`data' = `v'_`data'_95
			matrix `v'_`data' = `v'_`data', `v'_`data'_90[.,"ll"], `v'_`data'_90[.,"ul"] //* Adding 90% CI bounds from the 90% matrix (columns 4 and 5)

			* Rename columns for clarity
			matrix colnames `v'_`data' = b se t ll95 ul95 ll90 ul90

			* Verify the combined matrix
			matrix list `v'_`data'
			
			* Create AgSS matrix as a copy of LSMS structure
			matrix `v'_AgSS = `v'_FTF

			* Replace all entries in AgSS matrix with missing values
			local rows = rowsof(`v'_AgSS)
			local cols = colsof(`v'_AgSS)

			forvalues i = 1/`rows' {
				forvalues j = 1/`cols' {
					matrix `v'_AgSS[`i', `j'] = .
				}
			}
			
		}	
		
	}
	

	
	
	

*Figures - main manuscript
		
	//Figure 1: Timing of GoE fertilizer interventions, growing seasons, and surveys relative to the three analysis time periods. 
				*This figure is manually created.  

	//Figure 2: Effect of fertilizer blending facilities on fertilizer adoption: summary of main estimates by cohort and combined across cohorts.  
			*Fertilizer adpotion - binary outcome
				*This code preps the control group base period st dev to be added to the coeffplot
				cap drop _* S_n
				coefplot  ///
					(matrix(NPS_Dummy_FTF[,1]), label(NPS) ci((NPS_Dummy_FTF[.,4] NPS_Dummy_FTF[.,5] ) (NPS_Dummy_FTF[.,6] NPS_Dummy_FTF[.,7] )) )  ///
					(matrix(DAP_Dummy_FTF[,1]), label(DAP) ci((DAP_Dummy_FTF[.,4] DAP_Dummy_FTF[.,5] ) (DAP_Dummy_FTF[.,6] DAP_Dummy_FTF[.,7] )))  ///
					(matrix(Fert_Dummy_FTF[,1]), label(Overall fertilizer) ci((Fert_Dummy_FTF[.,4] Fert_Dummy_FTF[.,5] ) (Fert_Dummy_FTF[.,6] Fert_Dummy_FTF[.,7]))) ///
					,  generate 
						
					g __cgsdpos = .
					g __cgsdneg = .
					gen S_n=_n if __plot!=.
					loc i=1
					display 
					foreach fert in NPS DAP Fert {
						replace __cgsdpos = ``fert'_Dummy_max_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
						replace __cgsdneg = ``fert'_Dummy_min_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
						replace __cgsdpos = ``fert'_Dummy_max_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)
						replace __cgsdneg = ``fert'_Dummy_min_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)	
						loc i = `i'+1
						}	
				coefplot  ///
					(matrix(NPS_Dummy_FTF[,1]), label(NPS) ci((NPS_Dummy_FTF[.,4] NPS_Dummy_FTF[.,5] ) (NPS_Dummy_FTF[.,6] NPS_Dummy_FTF[.,7] )) )  ///
					(matrix(DAP_Dummy_FTF[,1]), label(DAP) ci((DAP_Dummy_FTF[.,4] DAP_Dummy_FTF[.,5] ) (DAP_Dummy_FTF[.,6] DAP_Dummy_FTF[.,7] )))  ///
					(matrix(Fert_Dummy_FTF[,1]), label(Overall fertilizer) ci((Fert_Dummy_FTF[.,4] Fert_Dummy_FTF[.,5] ) (Fert_Dummy_FTF[.,6] Fert_Dummy_FTF[.,7] ))) ///
					, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
					drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
					vertical legend(rows(1) size(small) region(lstyle(none)) pos(6))  order(3 "NPS" 6 "DAP" 9 "Mazie" 10 "Control group st dev")  ///
					graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
					msymbol(D)  levels( 95) ///
					yline(0, lstyle(1) lcolor(red)) yscale(r(-.75 .75)) ylabel(-.75(.25).75,nogrid) ///
					addplot((scatter __cgsdpos __at, msymbol(smplus)  mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus)  mcolor(gs8) msize(medium)))
					drop _* S_n
				graph export "Figures/Figure_2.png", replace

				
	//Figure 3: Effect of fertilizer blending facilities on fertilizer application rate (kg/ha): summary of main estimates by cohort and combined across cohorts.  
				coefplot  ///
					(matrix(w_NPS_rate_FTF[,1]), label("NPS") ci((w_NPS_rate_FTF[,4] w_NPS_rate_FTF[,5] ) (w_NPS_rate_FTF[,6] w_NPS_rate_FTF[,7] )) )  ///
					(matrix(w_DAP_rate_FTF[,1]), label("DAP") ci((w_DAP_rate_FTF[,4] w_DAP_rate_FTF[,5] ) (w_DAP_rate_FTF[,6] w_DAP_rate_FTF[,7] )))  ///
					(matrix(w_App_rate_FTF[,1]), label("Overall fertilizer") ci((w_App_rate_FTF[,4] w_App_rate_FTF[,5] ) (w_App_rate_FTF[,6] w_App_rate_FTF[,7]))) ///
					, generate replace

				g __cgsdpos = .
				g __cgsdneg = .
				gen S_n=_n if __plot!=.
				
				loc i=1
				display 
				foreach fert in NPS DAP App  {
					replace __cgsdpos = `w_`fert'_rate_max_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
					replace __cgsdneg = `w_`fert'_rate_min_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
					replace __cgsdpos = `w_`fert'_rate_max_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)
					replace __cgsdneg = `w_`fert'_rate_min_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)			
					loc i = `i'+1
					}
				coefplot  ///
					(matrix(w_NPS_rate_FTF[,1]), label("NPS") ci((w_NPS_rate_FTF[,4] w_NPS_rate_FTF[,5] ) (w_NPS_rate_FTF[,6] w_NPS_rate_FTF[,7] )))  ///
					(matrix(w_DAP_rate_FTF[,1]), label("DAP") ci((w_DAP_rate_FTF[,4] w_DAP_rate_FTF[,5] ) (w_DAP_rate_FTF[,6] w_DAP_rate_FTF[,7] )))  ///
					(matrix(w_App_rate_FTF[,1]), label("Overall fertilizer") ci((w_App_rate_FTF[,4] w_App_rate_FTF[,5] ) (w_App_rate_FTF[,6] w_App_rate_FTF[,7]))) ///
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(rows(1) size(small) region(lstyle(none)) order(3 "NPS" 6 "DAP" 9 "Overall fertilizer" 10 "Control group st dev"))    ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) legend(pos(6))  ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-50 50)) ylabel(-50(25)50,nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))		

				drop _* S_n
				graph export "Figures/Figure_3.png", replace

//Figure 4: Effect of fertilizer blending facilities on fertilizer (DAP) price: summary of main estimates by cohort and combined across cohorts.  
			
			coefplot  ///
				(matrix(DAP_Price_mr_FTF[,1]), label("") ci((DAP_Price_mr_FTF[.,4] DAP_Price_mr_FTF[.,5] ) (DAP_Price_mr_FTF[.,6] DAP_Price_mr_FTF[.,7] )) )  ///
				, generate replace
				
				g __cgsdpos = .
				g __cgsdneg = .
				gen S_n=_n if __plot!=.
				
				loc i=1
				display 
				foreach fert in DAP  {
					replace __cgsdpos = ``fert'_Price_mr_max_FTF_c1' if __plot==`i' & inrange(S_n, 1, 2)
					replace __cgsdneg = ``fert'_Price_mr_min_FTF_c1' if __plot==`i' & inrange(S_n, 1, 2)
					replace __cgsdpos = ``fert'_Price_mr_max_FTF_c2' if __plot==`i' & inlist(S_n, 3)
					replace __cgsdneg = ``fert'_Price_mr_min_FTF_c2' if __plot==`i' & inlist(S_n, 3)			
					loc i = `i'+1
					}
			coefplot  ///
				(matrix(DAP_Price_mr_FTF[,1]), label("") ci((DAP_Price_mr_FTF[.,4] DAP_Price_mr_FTF[.,5] ) (DAP_Price_mr_FTF[.,6] DAP_Price_mr_FTF[.,7] )) )  ///
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical  legend(off)   ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) legend(pos(6))  ///
				ytitle(Price of DAP in birr, size(medium)) ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-15 15)) ylabel(-15(5)15, nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))		
				drop _* S_n
				
				graph export "Figures/Figure_4.png", replace
			
			
	//Figure 5: Effect of fertilizer blending facilities on yield (kg/ha): summary of main estimates by cohort and combined across cohorts.  

				coefplot  ///
					(matrix(CL_w_Yield_kgha_1_FTF[,1]), label(Teff) ci((CL_w_Yield_kgha_1_FTF[.,4] CL_w_Yield_kgha_1_FTF[.,5] ) (CL_w_Yield_kgha_1_FTF[.,6] CL_w_Yield_kgha_1_FTF[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_2_FTF[,1]), label(wheat) ci((CL_w_Yield_kgha_2_FTF[.,4] CL_w_Yield_kgha_2_FTF[.,5] ) (CL_w_Yield_kgha_2_FTF[.,6] CL_w_Yield_kgha_2_FTF[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_3_FTF[,1]), label(Maize) ci((CL_w_Yield_kgha_3_FTF[.,4] CL_w_Yield_kgha_3_FTF[.,5] ) (CL_w_Yield_kgha_3_FTF[.,6] CL_w_Yield_kgha_3_FTF[.,7] ))) ///
					, generate replace
						
				g __cgsdpos = .
				g __cgsdneg = .
				gen S_n=_n if __plot!=.
				loc i=1
				display 
			foreach crop in 1 2 3  {
					replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
					replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_FTF_c1' if __plot==`i' & inrange(S_n, 1, 6)
					replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)
					replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_FTF_c2' if __plot==`i' & inrange(S_n, 7, 9)	
					loc i = `i'+1
					}	
				coefplot  ///
					(matrix(CL_w_Yield_kgha_1_FTF[,1]), label(Teff) ci((CL_w_Yield_kgha_1_FTF[.,4] CL_w_Yield_kgha_1_FTF[.,5] ) (CL_w_Yield_kgha_1_FTF[.,6] CL_w_Yield_kgha_1_FTF[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_2_FTF[,1]), label(wheat) ci((CL_w_Yield_kgha_2_FTF[.,4] CL_w_Yield_kgha_2_FTF[.,5] ) (CL_w_Yield_kgha_2_FTF[.,6] CL_w_Yield_kgha_2_FTF[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_3_FTF[,1]), label(Maize) ci((CL_w_Yield_kgha_3_FTF[.,4] CL_w_Yield_kgha_3_FTF[.,5] ) (CL_w_Yield_kgha_3_FTF[.,6] CL_w_Yield_kgha_3_FTF[.,7] ))) ///
						, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
						drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
						vertical legend(rows(1) size(small) region(lstyle(none)) pos(6))  order(3 "Teff" 6 "Wheat" 9 "Mazie" 10 "Control group st dev")  ///
						graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
						msymbol(D)  levels( 95) ///
						yline(0, lstyle(1) lcolor(red)) yscale(r(-1000 1000)) ylabel(-1000(500)1000,nogrid)	///	
						addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))
						drop _* S_n
						
						graph export "Figures/Figure_5.png", replace
				
						
	//Figure 6: Effect of fertilizer blending facilities on value of production (total and by crop) in birr: summary of main estimates by cohort and combined across cohorts.  
				coefplot  ///
					(matrix(TOT_V_harvest_kg_FTF[,1]), label(Total) ci((TOT_V_harvest_kg_FTF[.,4] TOT_V_harvest_kg_FTF[.,5] ) (TOT_V_harvest_kg_FTF[.,6] TOT_V_harvest_kg_FTF[.,7] )) )  ///
					(matrix(CL_V_w_harvest_kg_1_FTF[,1]), label(Teff) ci((CL_V_w_harvest_kg_1_FTF[.,4] CL_V_w_harvest_kg_1_FTF[.,5] ) (CL_V_w_harvest_kg_1_FTF[.,6] CL_V_w_harvest_kg_1_FTF[.,7] )) )  ///
					(matrix(CL_V_w_harvest_kg_2_FTF[,1]), label(wheat) ci((CL_V_w_harvest_kg_2_FTF[.,4] CL_V_w_harvest_kg_2_FTF[.,5] ) (CL_V_w_harvest_kg_2_FTF[.,6] CL_V_w_harvest_kg_2_FTF[.,7] )))  ///
					(matrix(CL_V_w_harvest_kg_3_FTF[,1]), label(Maize) ci((CL_V_w_harvest_kg_3_FTF[.,4] CL_V_w_harvest_kg_3_FTF[.,5] ) (CL_V_w_harvest_kg_3_FTF[.,6] CL_V_w_harvest_kg_3_FTF[.,7] ))) ///
					, generate replace

				
				g __cgsdpos = .
				g __cgsdneg = .
				gen S_n=_n if __plot!=.
				
				display 
					replace __cgsdpos = `TOT_V_harvest_kg_max_FTF_c1' if __plot==1 & inlist(S_n, 1, 5)
					replace __cgsdneg = `TOT_V_harvest_kg_min_FTF_c1' if __plot==1 & inlist(S_n, 1, 5)
					replace __cgsdpos = `TOT_V_harvest_kg_max_FTF_c2' if __plot==1 & inlist(S_n, 9)
					replace __cgsdneg = `TOT_V_harvest_kg_min_FTF_c2' if __plot==1 & inlist(S_n, 9)		
					loc i=2
				display 
				foreach crop in 1 2 3  {
					replace __cgsdpos = `CL_V_w_harvest_kg_`crop'_max_FTF_c1' if __plot==`i' & inlist(S_n, 2, 3, 4, 6, 7, 8)
					replace __cgsdneg = `CL_V_w_harvest_kg_`crop'_min_FTF_c1' if __plot==`i' & inlist(S_n,  2, 3, 4, 6, 7, 8)
					replace __cgsdpos = `CL_V_w_harvest_kg_`crop'_max_FTF_c2' if __plot==`i' & inrange(S_n, 10, 12)
					replace __cgsdneg = `CL_V_w_harvest_kg_`crop'_min_FTF_c2' if __plot==`i' & inrange(S_n, 10, 12)
					loc i = `i'+1
					}		
				coefplot (TOT_V_harvest_kg_FTF, label(Total) ) (CL_V_w_harvest_kg_1_FTF, label(Teff) ) (CL_V_w_harvest_kg_2_FTF, label(Wheat)) (CL_V_w_harvest_kg_3_FTF, label(Maize)) /// 
				, rename (GAverage= Combined G2 = "Cohort one" G3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(rows(1) size(small) region(lstyle(none)) pos(6))  order(3 "Total" 6 "Teff" 9 "Wheat" 12 "Mazie" 13 "Control group st dev")  ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) ///
				ytitle(Value of production in birr, size(medium)) ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-6000 6000)) ylabel(-6000(1500)6000,nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))		
				
				drop _* S_n				
				
				graph export "Figures/Figure_6.png", replace
						
			
		//Figure 7: Effect of fertilizer blending facilities on Consumption expenditure: summary of main estimates by cohort and combined across cohorts.  

			coefplot  ///
				(matrix(paexpD_R_FTF[,1]), label("Daily real per adult equivalent total expenditure (birr) in 2011 base price") ci((paexpD_R_FTF[.,4] paexpD_R_FTF[.,5] ) (paexpD_R_FTF[.,6] paexpD_R_FTF[.,7] )) )  ///
				, generate replace
				
				g __cgsdpos = .
				g __cgsdneg = .
				gen S_n=_n if __plot!=.
				
				loc i=1
				display 
				foreach exnd in paexpD_R  {
					replace __cgsdpos = ``exnd'_max_FTF_c1' if __plot==`i' & inrange(S_n, 1, 2)
					replace __cgsdneg = ``exnd'_min_FTF_c1' if __plot==`i' & inrange(S_n, 1, 2)
					replace __cgsdpos = ``exnd'_max_FTF_c2' if __plot==`i' & inlist(S_n, 3)
					replace __cgsdneg = ``exnd'_min_FTF_c2' if __plot==`i' & inlist(S_n, 3)			
					loc i = `i'+1
					}		
			coefplot  ///
				(matrix(paexpD_R_FTF[,1]), label("Daily real per adult equivalent total expenditure (birr) in 2011 base price") ci((paexpD_R_FTF[.,4] paexpD_R_FTF[.,5] ) (paexpD_R_FTF[.,6] paexpD_R_FTF[.,7] )) )  ///
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(off)  ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) ///
				ytitle(Daily real per adult equivalent total expenditure (birr) in 2011 base price, size(small)) ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-15 15)) ylabel(-15(5)15,nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))		
				
				drop _* S_n				
				
				graph export "Figures/Figure_7.png", replace	
				
	
	
	
	
	
*Figures - Supplemental Materials

	//Figure S.1: Average Ethiopia-wide fertilizer application rate, total fertilizer used, and number of extension package users in Ethiopia between 2005-06 and 2017-18.
				*Figure S.1. is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.
	//Figure S.2: 
			*(Panel a) Total Ethiopia-wide fertilizer use by product between 2013 and 2018
				*Figure S.2. Panle (a) is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.

			*(Panel b) share of fertilizer users using each product  between 2013 and 2018 - add
			preserve 
			use FTF.dta, clear 
				*Collapsing to get average (proportion) by year
				collapse (mean) DAP_Dummy Urea_Dummy NPS_Dummy, by(year)

				*Renaming for clarity 
				rename DAP_Dummy share_dap
				rename Urea_Dummy share_urea
				rename NPS_Dummy share_nps

				*Reshaping to long for easy plotting
				reshape long share_, i(year) j(fert_type) string
				rename share_ proportion

				*Ploting 
				twoway ///
					(connected proportion year if fert_type == "dap", msymbol(O) lcolor(blue) lwidth(medthick)) ///
					(connected proportion year if fert_type == "urea", msymbol(O) lcolor(orange) lwidth(medthick)) ///
					(connected proportion year if fert_type == "nps", msymbol(O) lcolor(gs10) lwidth(medthick)), ///
					legend(position(6) row(1) ///
						   label(1 "DAP") label(2 "Urea") label(3 "NPS")) ///
					ylabel(0(0.1)0.7) ///
					xlabel(1 "2013" 2 "2015" 3 "2018") ///
					ytitle("Proportions of fertilizer users") ///
					xtitle("Year") ///
					graphregion(color(white)) ///
					plotregion(color(white)) ///
					title("") ///
					name(fert_share_final, replace)

			restore 
			
	//Figure S.6: Effect of fertilizer blending facilities on fertilizer adoption: summary of estimates across datasets by cohort and combined across cohorts.  
				coefplot  ///
					(matrix(NPS_Dummy_FTF[,1]), label(NPS) ci((NPS_Dummy_FTF[.,4] NPS_Dummy_FTF[.,5] ) (NPS_Dummy_FTF[.,6] NPS_Dummy_FTF[.,7] )) )  ///
					(matrix(NPS_Dummy_AgSS[,1]), label(NPS) ci((NPS_Dummy_AgSS[.,4] NPS_Dummy_AgSS[.,5] ) (NPS_Dummy_AgSS[.,6] NPS_Dummy_AgSS[.,7] )) )  ///
					(matrix(NPS_Dummy_LSMS[,1]), label(NPS) ci((NPS_Dummy_LSMS[.,4] NPS_Dummy_LSMS[.,5] ) (NPS_Dummy_LSMS[.,6] NPS_Dummy_LSMS[.,7] )) )  ///
					(matrix(DAP_Dummy_FTF[,1]), label(DAP) ci((DAP_Dummy_FTF[.,4] DAP_Dummy_FTF[.,5] ) (DAP_Dummy_FTF[.,6] DAP_Dummy_FTF[.,7] )))  ///
					(matrix(DAP_Dummy_AgSS[,1]), label(DAP) ci((DAP_Dummy_AgSS[.,4] DAP_Dummy_AgSS[.,5] ) (DAP_Dummy_AgSS[.,6] DAP_Dummy_AgSS[.,7] )))  ///
					(matrix(DAP_Dummy_LSMS[,1]), label(DAP) ci((DAP_Dummy_LSMS[.,4] DAP_Dummy_LSMS[.,5] ) (DAP_Dummy_LSMS[.,6] DAP_Dummy_LSMS[.,7] )))  ///
					(matrix(Fert_Dummy_FTF[,1]), label(Overall fertilizer) ci((Fert_Dummy_FTF[.,4] Fert_Dummy_FTF[.,5] ) (Fert_Dummy_FTF[.,6] Fert_Dummy_FTF[.,7]))) ///
					(matrix(Fert_Dummy_AgSS[,1]), label(Overall fertilizer) ci((Fert_Dummy_AgSS[.,4] Fert_Dummy_AgSS[.,5] ) (Fert_Dummy_AgSS[.,6] Fert_Dummy_AgSS[.,7]))) ///
					(matrix(Fert_Dummy_LSMS[,1]), label(Overall fertilizer) ci((Fert_Dummy_LSMS[.,4] Fert_Dummy_LSMS[.,5] ) (Fert_Dummy_LSMS[.,6] Fert_Dummy_LSMS[.,7]))) ///
					, generate replace
					
			gen index = cond(!missing(__plot), _n, .)

			* Generate placeholders for control group standard deviations
			gen __cgsdpos = .
			gen __cgsdneg = .

			* Initialize a counter for the plot index 
			local i = 1
			* Loop through fertilizer types and datasets to fill the scatterplot values according to coding above
			foreach fert in NPS DAP Fert {
				foreach data in FTF AgSS LSMS {
					replace __cgsdpos = ``fert'_Dummy_max_`data'_c1' if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdneg = ``fert'_Dummy_min_`data'_c1' if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdpos = ``fert'_Dummy_max_`data'_c2' if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdneg = ``fert'_Dummy_min_`data'_c2' if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdpos = . if  inlist(index, 3, 6, 9)
					replace __cgsdneg = . if  inlist(index, 3, 6, 9)
					loc i = `i'+1
					}
				}

				
			* Code to manually overwrite issue with addplot [set s start value in line 113 to where things are cutting off]
			count if __at != . 
			loc S = r(N)
			g __at2 = .
			g __cgsdpos2 = .
			g __cgsdneg2 = .
			loc i=1
			forvalues s=19/`S' {
				su __at if _n==`s'
				loc val = r(mean)
				su __cgsdpos if _n==`s'
				loc valpos = r(mean)
				su __cgsdneg if _n==`s'
				loc valneg = r(mean)
				replace __at2 = `val' in `i'
				replace __cgsdpos2 = `valpos' in `i'
				replace __cgsdneg2 = `valneg' in `i'
				loc i = `i'+1
				}
			
			coefplot  ///
					(matrix(NPS_Dummy_FTF[,1]), label("") msymbol(O) pstyle(p1) ci((NPS_Dummy_FTF[.,4] NPS_Dummy_FTF[.,5] ) (NPS_Dummy_FTF[.,6] NPS_Dummy_FTF[.,7] )) )  ///
					(matrix(NPS_Dummy_AgSS[,1]), label("") msymbol(S) pstyle(p1) ci((NPS_Dummy_AgSS[.,4] NPS_Dummy_AgSS[.,5] ) (NPS_Dummy_AgSS[.,6] NPS_Dummy_AgSS[.,7] )) )  ///
					(matrix(NPS_Dummy_LSMS[,1]), label("") msymbol(T) pstyle(p1)  drop (r1) ci((NPS_Dummy_LSMS[.,4] NPS_Dummy_LSMS[.,5] ) (NPS_Dummy_LSMS[.,6] NPS_Dummy_LSMS[.,7] )) )  ///
					(matrix(DAP_Dummy_FTF[,1]), label("") msymbol(O) pstyle(p2) ci((DAP_Dummy_FTF[.,4] DAP_Dummy_FTF[.,5] ) (DAP_Dummy_FTF[.,6] DAP_Dummy_FTF[.,7] )))  ///
					(matrix(DAP_Dummy_AgSS[,1]), label("") msymbol(S) pstyle(p2) ci((DAP_Dummy_AgSS[.,4] DAP_Dummy_AgSS[.,5] ) (DAP_Dummy_AgSS[.,6] DAP_Dummy_AgSS[.,7] )))  ///
					(matrix(DAP_Dummy_LSMS[,1]), label("") msymbol(T) pstyle(p2) drop (r1) ci((DAP_Dummy_LSMS[.,4] DAP_Dummy_LSMS[.,5] ) (DAP_Dummy_LSMS[.,6] DAP_Dummy_LSMS[.,7] )))  ///
					(matrix(Fert_Dummy_FTF[,1]), label("") msymbol(O) pstyle(p3) ci((Fert_Dummy_FTF[.,4] Fert_Dummy_FTF[.,5] ) (Fert_Dummy_FTF[.,6] Fert_Dummy_FTF[.,7]))) ///
					(matrix(Fert_Dummy_AgSS[,1]), label("") msymbol(S) pstyle(p3) ci((Fert_Dummy_AgSS[.,4] Fert_Dummy_AgSS[.,5] ) (Fert_Dummy_AgSS[.,6] Fert_Dummy_AgSS[.,7]))) ///
					(matrix(Fert_Dummy_LSMS[,1]), label("") msymbol(T) pstyle(p3) drop (r1) ci((Fert_Dummy_LSMS[.,4] Fert_Dummy_LSMS[.,5] ) (Fert_Dummy_LSMS[.,6] Fert_Dummy_LSMS[.,7]))) ///
			,   rename (r1= Combined r2 = "Cohort one" r3="Cohort two") ///
				vertical legend(rows(5) pos(6) size(small) region(lstyle(none)) /// 
				order(8 "NPS"   17 "DAP" 20 "Overall fertilizer"  3 12 21 "FTF"   6 15 24 "AgSS- not displayed" 9 18 27 "LSMS" 29 "" 29 "" 29 "Control group st dev")) /// 
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients") ///
				levels(95 90) yline(0, lstyle(1) lcolor(red)) ///
				yscale(r(-.75 .75)) ylabel(-.75(.25).75, nogrid) ///
				addplot((scatter __cgsdpos __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdpos2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)))

				drop _* index
				
			* Export the combined graph	
				
				graph export "Figures/Figure_S6.png", replace
				

	// Figure S.7: Effect of fertilizer blending facilities on fertilizer application rate (kg/ha): summary of estimates across datasets by cohort and combined across cohorts.  
			coefplot  ///
					(matrix(w_NPS_rate_FTF[,1]),  ci((w_NPS_rate_FTF[.,4] w_NPS_rate_FTF[.,5] ) (w_NPS_rate_FTF[.,6] w_NPS_rate_FTF[.,7] )) )  ///
					(matrix(w_NPS_rate_AgSS[,1]), ci((w_NPS_rate_AgSS[.,4] w_NPS_rate_AgSS[.,5] ) (w_NPS_rate_AgSS[.,6] w_NPS_rate_AgSS[.,7] )) )  ///
					(matrix(w_NPS_rate_LSMS[,1]), ci((w_NPS_rate_LSMS[.,4] w_NPS_rate_LSMS[.,5] ) (w_NPS_rate_LSMS[.,6] w_NPS_rate_LSMS[.,7] )) )  ///
					(matrix(w_DAP_rate_FTF[,1]),  ci((w_DAP_rate_FTF[.,4] w_DAP_rate_FTF[.,5] ) (w_DAP_rate_FTF[.,6] w_DAP_rate_FTF[.,7] )))  ///
					(matrix(w_DAP_rate_AgSS[,1]), ci((w_DAP_rate_AgSS[.,4] w_DAP_rate_AgSS[.,5] ) (w_DAP_rate_AgSS[.,6] w_DAP_rate_AgSS[.,7] )))  ///
					(matrix(w_DAP_rate_LSMS[,1]), ci((w_DAP_rate_LSMS[.,4] w_DAP_rate_LSMS[.,5] ) (w_DAP_rate_LSMS[.,6] w_DAP_rate_LSMS[.,7] )))  ///
					(matrix(w_App_rate_FTF[,1]),  ci((w_App_rate_FTF[.,4] w_App_rate_FTF[.,5] ) (w_App_rate_FTF[.,6] w_App_rate_FTF[.,7]))) ///
					(matrix(w_App_rate_AgSS[,1]), ci((w_App_rate_AgSS[.,4] w_App_rate_AgSS[.,5] ) (w_App_rate_AgSS[.,6] w_App_rate_AgSS[.,7]))) ///
					(matrix(w_App_rate_LSMS[,1]), ci((w_App_rate_LSMS[.,4] w_App_rate_LSMS[.,5] ) (w_App_rate_LSMS[.,6] w_App_rate_LSMS[.,7]))) ///
					, generate replace
			gen index = cond(!missing(__plot), _n, .)

			* Generate placeholders for control group standard deviations
			gen __cgsdpos = .
			gen __cgsdneg = .

			* Initialize a counter for the plot index 
			local i = 1
			* Loop through fertilizer types and datasets to fill the scatterplot values according to coding above
			foreach fert in NPS DAP App {
				foreach data in FTF AgSS LSMS {
					replace __cgsdpos = `w_`fert'_rate_max_`data'_c1' if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdneg = `w_`fert'_rate_min_`data'_c1' if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdpos = `w_`fert'_rate_max_`data'_c2' if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdneg = `w_`fert'_rate_min_`data'_c2' if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdpos = . if  inlist(index, 3, 6, 9)
					replace __cgsdneg = . if  inlist(index, 3, 6, 9)
					loc i = `i'+1
					}
				}
				
			** Code to manually overwrite issue with addplot [set s start value in line 113 to where things are cutting off]
			count if __at != . 
			loc S = r(N)
			g __at2 = .
			g __cgsdpos2 = .
			g __cgsdneg2 = .
			loc i=1
			forvalues s=24/`S' {
				su __at if _n==`s'
				loc val = r(mean)
				su __cgsdpos if _n==`s'
				loc valpos = r(mean)
				su __cgsdneg if _n==`s'
				loc valneg = r(mean)
				replace __at2 = `val' in `i'
				replace __cgsdpos2 = `valpos' in `i'
				replace __cgsdneg2 = `valneg' in `i'
				loc i = `i'+1
				}

			coefplot  ///
					(matrix(w_NPS_rate_FTF[,1]), label("") msymbol(O) pstyle(p1) ci((w_NPS_rate_FTF[.,4] w_NPS_rate_FTF[.,5] ) (w_NPS_rate_FTF[.,6] w_NPS_rate_FTF[.,7] )) )  ///
					(matrix(w_NPS_rate_AgSS[,1]), label("") msymbol(S) pstyle(p1) ci((w_NPS_rate_AgSS[.,4] w_NPS_rate_AgSS[.,5] ) (w_NPS_rate_AgSS[.,6] w_NPS_rate_AgSS[.,7] )) )  ///
					(matrix(w_NPS_rate_LSMS[,1]), label("") msymbol(T) pstyle(p1)  drop (r1) ci((w_NPS_rate_LSMS[.,4] w_NPS_rate_LSMS[.,5] ) (w_NPS_rate_LSMS[.,6] w_NPS_rate_LSMS[.,7] )) )  ///
					(matrix(w_DAP_rate_FTF[,1]), label("") msymbol(O) pstyle(p2) ci((w_DAP_rate_FTF[.,4] w_DAP_rate_FTF[.,5] ) (w_DAP_rate_FTF[.,6] w_DAP_rate_FTF[.,7] )))  ///
					(matrix(w_DAP_rate_AgSS[,1]), label("") msymbol(S) pstyle(p2) ci((w_DAP_rate_AgSS[.,4] w_DAP_rate_AgSS[.,5] ) (w_DAP_rate_AgSS[.,6] w_DAP_rate_AgSS[.,7] )))  ///
					(matrix(w_DAP_rate_LSMS[,1]), label("") msymbol(T) pstyle(p2) drop (r1) ci((w_DAP_rate_LSMS[.,4] w_DAP_rate_LSMS[.,5] ) (w_DAP_rate_LSMS[.,6] w_DAP_rate_LSMS[.,7] )))  ///
					(matrix(w_App_rate_FTF[,1]), label("") msymbol(O) pstyle(p3) ci((w_App_rate_FTF[.,4] w_App_rate_FTF[.,5] ) (w_App_rate_FTF[.,6] w_App_rate_FTF[.,7]))) ///
					(matrix(w_App_rate_AgSS[,1]), label("") msymbol(S) pstyle(p3) ci((w_App_rate_AgSS[.,4] w_App_rate_AgSS[.,5] ) (w_App_rate_AgSS[.,6] w_App_rate_AgSS[.,7]))) ///
					(matrix(w_App_rate_LSMS[,1]), label("") msymbol(T) pstyle(p3) drop (r1) ci((w_App_rate_LSMS[.,4] w_App_rate_LSMS[.,5] ) (w_App_rate_LSMS[.,6] w_App_rate_LSMS[.,7]))) ///
			,   rename (r1= Combined r2 = "Cohort one" r3="Cohort two") ///
				vertical legend(rows(5) pos(6) size(small) region(lstyle(none)) /// 
				order(8 "NPS"   17 "DAP" 20 "Overall fertilizer"  3 12 21 "FTF"   6 15 24 "AgSS- not displayed"" 9 18 27 "LSMS" 29 "" 29 "" 29 "Control group st dev")) /// 
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients") ///
				levels(95 90) yline(0, lstyle(1) lcolor(red)) ///
				yscale(r(-250 250)) ylabel(-250(100)250, nogrid) ///
				addplot((scatter __cgsdpos __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdpos2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)))

			drop _* index
			
			* Export the combined graph

			graph export "Figures/Figure_S7.png", replace
			
	// Figure S.8: Effect of fertilizer blending facilities on yield (kg/ha): summary of estimates across datasets by cohort and combined across cohorts.  
			coefplot  ///
					(matrix(CL_w_Yield_kgha_1_FTF[,1]),  ci((CL_w_Yield_kgha_1_FTF[.,4] CL_w_Yield_kgha_1_FTF[.,5] ) (CL_w_Yield_kgha_1_FTF[.,6] CL_w_Yield_kgha_1_FTF[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_1_AgSS[,1]), ci((CL_w_Yield_kgha_1_AgSS[.,4] CL_w_Yield_kgha_1_AgSS[.,5] ) (CL_w_Yield_kgha_1_AgSS[.,6] CL_w_Yield_kgha_1_AgSS[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_1_LSMS[,1]), ci((CL_w_Yield_kgha_1_LSMS[.,4] CL_w_Yield_kgha_1_LSMS[.,5] ) (CL_w_Yield_kgha_1_LSMS[.,6] CL_w_Yield_kgha_1_LSMS[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_2_FTF[,1]),  ci((CL_w_Yield_kgha_2_FTF[.,4] CL_w_Yield_kgha_2_FTF[.,5] ) (CL_w_Yield_kgha_2_FTF[.,6] CL_w_Yield_kgha_2_FTF[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_2_AgSS[,1]), ci((CL_w_Yield_kgha_2_AgSS[.,4] CL_w_Yield_kgha_2_AgSS[.,5] ) (CL_w_Yield_kgha_2_AgSS[.,6] CL_w_Yield_kgha_2_AgSS[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_2_LSMS[,1]), ci((CL_w_Yield_kgha_2_LSMS[.,4] CL_w_Yield_kgha_2_LSMS[.,5] ) (CL_w_Yield_kgha_2_LSMS[.,6] CL_w_Yield_kgha_2_LSMS[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_3_FTF[,1]),  ci((CL_w_Yield_kgha_3_FTF[.,4] CL_w_Yield_kgha_3_FTF[.,5] ) (CL_w_Yield_kgha_3_FTF[.,6] CL_w_Yield_kgha_3_FTF[.,7]))) ///
					(matrix(CL_w_Yield_kgha_3_AgSS[,1]), ci((CL_w_Yield_kgha_3_AgSS[.,4] CL_w_Yield_kgha_3_AgSS[.,5] ) (CL_w_Yield_kgha_3_AgSS[.,6] CL_w_Yield_kgha_3_AgSS[.,7]))) ///
					(matrix(CL_w_Yield_kgha_3_LSMS[,1]), ci((CL_w_Yield_kgha_3_LSMS[.,4] CL_w_Yield_kgha_3_LSMS[.,5] ) (CL_w_Yield_kgha_3_LSMS[.,6] CL_w_Yield_kgha_3_LSMS[.,7]))) ///
					, generate replace
			gen index = cond(!missing(__plot), _n, .)

			* Generate placeholders for control group standard deviations
			gen __cgsdpos = .
			gen __cgsdneg = .

			* Initialize a counter for the plot index 
			local i = 1
			* Loop through fertilizer types and datasets to fill the scatterplot values according to coding above
			foreach crop in 1 2 3 {
				foreach data in FTF AgSS LSMS {
					replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_`data'_c1'  if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_`data'_c1' if __plot == `i' & inrange(index, 1, 18)
					replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_`data'_c2'  if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_`data'_c2' if __plot == `i' & inrange(index, 19, 24)
					replace __cgsdpos = . if  inlist(index, 3, 6, 9)
					replace __cgsdneg = . if  inlist(index, 3, 6, 9)
					loc i = `i'+1
					}
				}
				
			** Code to manually overwrite issue with addplot [set s start value in line 113 to where things are cutting off]
			count if __at != . 
			loc S = r(N)
			g __at2 = .
			g __cgsdpos2 = .
			g __cgsdneg2 = .
			loc i=1
			forvalues s=19/`S' {
				su __at if _n==`s'
				loc val = r(mean)
				su __cgsdpos if _n==`s'
				loc valpos = r(mean)
				su __cgsdneg if _n==`s'
				loc valneg = r(mean)
				replace __at2 = `val' in `i'
				replace __cgsdpos2 = `valpos' in `i'
				replace __cgsdneg2 = `valneg' in `i'
				loc i = `i'+1
				}

			coefplot  ///
					(matrix(CL_w_Yield_kgha_1_FTF[,1]), label("") msymbol(O) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF[.,4] CL_w_Yield_kgha_1_FTF[.,5] ) (CL_w_Yield_kgha_1_FTF[.,6] CL_w_Yield_kgha_1_FTF[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_1_AgSS[,1]), label("") msymbol(S) pstyle(p1) ci((CL_w_Yield_kgha_1_AgSS[.,4] CL_w_Yield_kgha_1_AgSS[.,5] ) (CL_w_Yield_kgha_1_AgSS[.,6] CL_w_Yield_kgha_1_AgSS[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_1_LSMS[,1]), label("") msymbol(T) pstyle(p1)  drop (r1) ci((CL_w_Yield_kgha_1_LSMS[.,4] CL_w_Yield_kgha_1_LSMS[.,5] ) (CL_w_Yield_kgha_1_LSMS[.,6] CL_w_Yield_kgha_1_LSMS[.,7] )) )  ///
					(matrix(CL_w_Yield_kgha_2_FTF[,1]), label("") msymbol(O) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF[.,4] CL_w_Yield_kgha_2_FTF[.,5] ) (CL_w_Yield_kgha_2_FTF[.,6] CL_w_Yield_kgha_2_FTF[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_2_AgSS[,1]), label("") msymbol(S) pstyle(p2) ci((CL_w_Yield_kgha_2_AgSS[.,4] CL_w_Yield_kgha_2_AgSS[.,5] ) (CL_w_Yield_kgha_2_AgSS[.,6] CL_w_Yield_kgha_2_AgSS[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_2_LSMS[,1]), label("") msymbol(T) pstyle(p2) drop (r1) ci((CL_w_Yield_kgha_2_LSMS[.,4] CL_w_Yield_kgha_2_LSMS[.,5] ) (CL_w_Yield_kgha_2_LSMS[.,6] CL_w_Yield_kgha_2_LSMS[.,7] )))  ///
					(matrix(CL_w_Yield_kgha_3_FTF[,1]), label("") msymbol(O) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF[.,4] CL_w_Yield_kgha_3_FTF[.,5] ) (CL_w_Yield_kgha_3_FTF[.,6] CL_w_Yield_kgha_3_FTF[.,7]))) ///
					(matrix(CL_w_Yield_kgha_3_AgSS[,1]), label("") msymbol(S) pstyle(p3) ci((CL_w_Yield_kgha_3_AgSS[.,4] CL_w_Yield_kgha_3_AgSS[.,5] ) (CL_w_Yield_kgha_3_AgSS[.,6] CL_w_Yield_kgha_3_AgSS[.,7]))) ///
					(matrix(CL_w_Yield_kgha_3_LSMS[,1]), label("") msymbol(T) pstyle(p3) drop (r1) ci((CL_w_Yield_kgha_3_LSMS[.,4] CL_w_Yield_kgha_3_LSMS[.,5] ) (CL_w_Yield_kgha_3_LSMS[.,6] CL_w_Yield_kgha_3_LSMS[.,7]))) ///
			,   rename (r1= Combined r2 = "Cohort one" r3="Cohort two") ///
				vertical legend(rows(5) pos(6) size(small) region(lstyle(none)) /// 
				order(8 "Teff"   17 "Wheat" 20 "Maize"  3 12 21 "FTF"   6 15 24 "AgSS- not displayed"" 9 18 27 "LSMS" 29 "" 29 "" 29 "Control group st dev")) /// 
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients") ///
				levels(95 90) yline(0, lstyle(1) lcolor(red)) ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-1200 1200)) ylabel(-1200(400)1200,nogrid)	///	
				addplot((scatter __cgsdpos __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg __at , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdpos2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)) ///
						(scatter __cgsdneg2 __at2 , msymbol(smplus) mcolor(gs8) msize(medium)))

			drop _* index
			
			* Export the combined graph
			graph export "Figures/Figure_S8.png", replace
			*/
	//Figure S.9: Pre-trend tests for the DAP and overall fertilizer adoption outcomes.
				*Figure S.9 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.
	//Figure S.10: Pre-trend tests for the DAP and overall fertilizer application rate outcomes.
				*Figure S.10 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.
	//Figure S.11: Pre-trend tests for the teff, wheat and maize yield outcomes.
				*Figure S.11 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.
		//Figure S.12: Effect of FBFs on fertilizer adoption: summary of estimates using different travel times to define FBF service areas.

	//Figure S.12, S.13, and S.14
			*Robustness - different travel times to define FBF service areas - estimations. 
			use FTF.dta, clear 
			keep if inlist(real(substr(string(A05, "%9.0f"), 1, 1)), 1, 3, 4, 7) //Restricting AgSS and LSMS datasets to the regiosn of the FTF data 
			
			*Defining outcomes for robustness - different cut off points
			loc outcome_robustness NPS_Dummy DAP_Dummy Fert_Dummy /// // Fertilizer adoption - binary outcome   
			w_App_rate w_DAP_rate w_NPS_rate /// //Application rates - all plots 
			CL_w_Yield_kgha_1 CL_w_Yield_kgha_2 CL_w_Yield_kgha_3  //Yield - crop level	
			

			*Define covariates 
			loc covariates land_area_sqm_10thou Female Litrate  
			foreach v in `outcome_robustness' {
				* Run csdid for each definition of treatment with for 90 and 95 % CI separately 
				foreach d in 5 4 3 2 {
				csdid `v' `covariates', ivar(id) time(year) gvar(first_treat`d') method(dripw) notyet wboot agg(group) rseed(19810310) cluster(A05) level(95)
				
				*Storing control group sds by cohort			
				estimates store `v'_FTF_`d'
				
				su `v' if (first_treat`d'==0| first_treat`d'==3) & year==1
			
				loc `v'_max_FTF_`d'_c1 = r(sd)
				display "`v'_max_FTF_`d'_c1: ``v'_max_FTF_`d'_c1"
				loc `v'_min_FTF_`d'_c1 = ``v'_max_FTF_`d'_c1'*-1
				
				su `v' if first_treat`d'==0 & year==2
				loc `v'_max_FTF_`d'_c2 = r(sd)
				display "`v'_max_FTF_c2: ``v'_max_FTF_c2'"
				loc `v'_min_FTF_`d'_c2 = ``v'_max_FTF_`d'_c2'*-1			
				
				*Storing 95% CIs 				
				matrix `v'_`d'_FTF_95 = e(cband)
				
				csdid `v' `covariates', ivar(id) time(year) gvar(first_treat`d') method(dripw) notyet wboot agg(group) rseed(19810310) cluster(A05) level(90)
					
				*Storing 90% CIs 				
				matrix `v'_`d'_FTF_90 = e(cband)
						
				*Combining the 90 and 95 % CIs. 
				matrix `v'_FTF_`d' = `v'_`d'_FTF_95
				matrix `v'_FTF_`d' = `v'_FTF_`d', `v'_`d'_FTF_90[.,"ll"], `v'_`d'_FTF_90[.,"ul"] //* Adding 90% CI bounds from the 90% matrix (columns 4 and 5)

				* Rename columns for clarity
				matrix colnames `v'_FTF_`d' = b se t ll95 ul95 ll90 ul90

				* Verify the combined matrix
				matrix list `v'_FTF_`d'
			}	
			}
		
	*Ploting results 
	//Figure S.12: Effect of fertilizer blending facilities on fertilizer adoption: summary of estimates across different cut off points for treatment by cohort and combined across cohorts.  	*Fertilizer adpotion - binary outcome
				*This code preps the control group base period st dev to be added to the coeffplot
				coefplot  ///
				(matrix(NPS_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((NPS_Dummy_FTF_5[.,4] NPS_Dummy_FTF_5[.,5] ) (NPS_Dummy_FTF_5[.,6] NPS_Dummy_FTF_5[.,7] ))) ///
				(matrix(NPS_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((NPS_Dummy_FTF_4[.,4] NPS_Dummy_FTF_4[.,5] ) (NPS_Dummy_FTF_4[.,6] NPS_Dummy_FTF_4[.,7] ))) ///		
				(matrix(NPS_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((NPS_Dummy_FTF_3[.,4] NPS_Dummy_FTF_3[.,5] ) (NPS_Dummy_FTF_3[.,6] NPS_Dummy_FTF_3[.,7] ))) ///	
				(matrix(NPS_Dummy_FTF_2[,1]), label("") msymbol(X) pstyle(p1) ci((NPS_Dummy_FTF_2[.,4] NPS_Dummy_FTF_2[.,5] ) (NPS_Dummy_FTF_2[.,6] NPS_Dummy_FTF_2[.,7] ))) ///
				(matrix(DAP_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((DAP_Dummy_FTF_5[.,4] DAP_Dummy_FTF_5[.,5] ) (DAP_Dummy_FTF_5[.,6] DAP_Dummy_FTF_5[.,7] ))) ///
				(matrix(DAP_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((DAP_Dummy_FTF_4[.,4] DAP_Dummy_FTF_4[.,5] ) (DAP_Dummy_FTF_4[.,6] DAP_Dummy_FTF_4[.,7] ))) ///		
				(matrix(DAP_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((DAP_Dummy_FTF_3[.,4] DAP_Dummy_FTF_3[.,5] ) (DAP_Dummy_FTF_3[.,6] DAP_Dummy_FTF_3[.,7] ))) ///	
				(matrix(DAP_Dummy_FTF_2[,1]), label("") msymbol(X) pstyle(p1) ci((DAP_Dummy_FTF_2[.,4] DAP_Dummy_FTF_2[.,5] ) (DAP_Dummy_FTF_2[.,6] DAP_Dummy_FTF_2[.,7] ))) ///	
				(matrix(Fert_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((Fert_Dummy_FTF_5[.,4] Fert_Dummy_FTF_5[.,5] ) (Fert_Dummy_FTF_5[.,6] Fert_Dummy_FTF_5[.,7] ))) ///
				(matrix(Fert_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((Fert_Dummy_FTF_4[.,4] Fert_Dummy_FTF_4[.,5] ) (Fert_Dummy_FTF_4[.,6] Fert_Dummy_FTF_4[.,7] ))) ///		
				(matrix(Fert_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((Fert_Dummy_FTF_3[.,4] Fert_Dummy_FTF_3[.,5] ) (Fert_Dummy_FTF_3[.,6] Fert_Dummy_FTF_3[.,7] ))) ///	
				(matrix(Fert_Dummy_FTF_2[,1]), label("") msymbol(X) pstyle(p1) ci((Fert_Dummy_FTF_2[.,4] Fert_Dummy_FTF_2[.,5] ) (Fert_Dummy_FTF_2[.,6] Fert_Dummy_FTF_2[.,7] ))) ///		
					, generate replace
				
				g __cgsdpos = .
				g __cgsdneg = .
				gen index = cond(!missing(__plot), _n, .)
				
				loc i=1
				display 
				foreach fert in NPS DAP Fert {
					foreach d in 5 4 3 2 {
					replace __cgsdpos = ``fert'_Dummy_max_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
					replace __cgsdneg = ``fert'_Dummy_min_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
					replace __cgsdpos = ``fert'_Dummy_max_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
					replace __cgsdneg = ``fert'_Dummy_min_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
					loc i = `i'+1
					}
					}
				coefplot  ///
				(matrix(NPS_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((NPS_Dummy_FTF_5[.,4] NPS_Dummy_FTF_5[.,5] ) (NPS_Dummy_FTF_5[.,6] NPS_Dummy_FTF_5[.,7] ))) ///
				(matrix(NPS_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((NPS_Dummy_FTF_4[.,4] NPS_Dummy_FTF_4[.,5] ) (NPS_Dummy_FTF_4[.,6] NPS_Dummy_FTF_4[.,7] ))) ///		
				(matrix(NPS_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((NPS_Dummy_FTF_3[.,4] NPS_Dummy_FTF_3[.,5] ) (NPS_Dummy_FTF_3[.,6] NPS_Dummy_FTF_3[.,7] ))) ///	
				(matrix(NPS_Dummy_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p1) ci((NPS_Dummy_FTF_2[.,4] NPS_Dummy_FTF_2[.,5] ) (NPS_Dummy_FTF_2[.,6] NPS_Dummy_FTF_2[.,7] ))) ///
				(matrix(DAP_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p2) ci((DAP_Dummy_FTF_5[.,4] DAP_Dummy_FTF_5[.,5] ) (DAP_Dummy_FTF_5[.,6] DAP_Dummy_FTF_5[.,7] ))) ///
				(matrix(DAP_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p2) ci((DAP_Dummy_FTF_4[.,4] DAP_Dummy_FTF_4[.,5] ) (DAP_Dummy_FTF_4[.,6] DAP_Dummy_FTF_4[.,7] ))) ///		
				(matrix(DAP_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p2) ci((DAP_Dummy_FTF_3[.,4] DAP_Dummy_FTF_3[.,5] ) (DAP_Dummy_FTF_3[.,6] DAP_Dummy_FTF_3[.,7] ))) ///	
				(matrix(DAP_Dummy_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p2) ci((DAP_Dummy_FTF_2[.,4] DAP_Dummy_FTF_2[.,5] ) (DAP_Dummy_FTF_2[.,6] DAP_Dummy_FTF_2[.,7] ))) ///	
				(matrix(Fert_Dummy_FTF_5[,1]), label("") msymbol(O) pstyle(p3) ci((Fert_Dummy_FTF_5[.,4] Fert_Dummy_FTF_5[.,5] ) (Fert_Dummy_FTF_5[.,6] Fert_Dummy_FTF_5[.,7] ))) ///
				(matrix(Fert_Dummy_FTF_4[,1]), label("") msymbol(S) pstyle(p3) ci((Fert_Dummy_FTF_4[.,4] Fert_Dummy_FTF_4[.,5] ) (Fert_Dummy_FTF_4[.,6] Fert_Dummy_FTF_4[.,7] ))) ///		
				(matrix(Fert_Dummy_FTF_3[,1]), label("") msymbol(T) pstyle(p3) ci((Fert_Dummy_FTF_3[.,4] Fert_Dummy_FTF_3[.,5] ) (Fert_Dummy_FTF_3[.,6] Fert_Dummy_FTF_3[.,7] ))) ///	
				(matrix(Fert_Dummy_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p3) ci((Fert_Dummy_FTF_2[.,4] Fert_Dummy_FTF_2[.,5] ) (Fert_Dummy_FTF_2[.,6] Fert_Dummy_FTF_2[.,7] ))) ///	
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(rows(6) pos(6) size(small) region(lstyle(none)) order(2 "NPS" 14 "DAP" 26 "Overall Fertilizer"  3 15 27 ///
				"5 hours" 6 18 30 "4 hours" 9 21 33 "3 hours" 12 24 36 "2 hours"  37 "" 37 "" 37 "Control group st dev"))    ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) legend(pos(6))  ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-.75 .75)) ylabel(-.75(.25).75,nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))
				
				drop _* index
				
				graph export "Figures/Figure_S12.png", replace
				
	//Figure S.13: Effect of fertilizer blending facilities on fertilizer application rate (kg/ha): summary of estimates across different cut off points for treatment by cohort and combined across cohorts.  	*Fertilizer adpotion - binary outcome
				*This code preps the control group base period st dev to be added to the coeffplot
			coefplot  ///
				(matrix(w_NPS_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((w_NPS_rate_FTF_5[.,4] w_NPS_rate_FTF_5[.,5] ) (w_NPS_rate_FTF_5[.,6] w_NPS_rate_FTF_5[.,7] ))) ///
				(matrix(w_NPS_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((w_NPS_rate_FTF_4[.,4] w_NPS_rate_FTF_4[.,5] ) (w_NPS_rate_FTF_4[.,6] w_NPS_rate_FTF_4[.,7] ))) ///		
				(matrix(w_NPS_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((w_NPS_rate_FTF_3[.,4] w_NPS_rate_FTF_3[.,5] ) (w_NPS_rate_FTF_3[.,6] w_NPS_rate_FTF_3[.,7] ))) ///	
				(matrix(w_NPS_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p1) ci((w_NPS_rate_FTF_2[.,4] w_NPS_rate_FTF_2[.,5] ) (w_NPS_rate_FTF_2[.,6] w_NPS_rate_FTF_2[.,7] ))) ///
				(matrix(w_DAP_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p2) ci((w_DAP_rate_FTF_5[.,4] w_DAP_rate_FTF_5[.,5] ) (w_DAP_rate_FTF_5[.,6] w_DAP_rate_FTF_5[.,7] ))) ///
				(matrix(w_DAP_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p2) ci((w_DAP_rate_FTF_4[.,4] w_DAP_rate_FTF_4[.,5] ) (w_DAP_rate_FTF_4[.,6] w_DAP_rate_FTF_4[.,7] ))) ///		
				(matrix(w_DAP_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p2) ci((w_DAP_rate_FTF_3[.,4] w_DAP_rate_FTF_3[.,5] ) (w_DAP_rate_FTF_3[.,6] w_DAP_rate_FTF_3[.,7] ))) ///	
				(matrix(w_DAP_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p2) ci((w_DAP_rate_FTF_2[.,4] w_DAP_rate_FTF_2[.,5] ) (w_DAP_rate_FTF_2[.,6] w_DAP_rate_FTF_2[.,7] ))) ///	
				(matrix(w_App_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p3) ci((w_App_rate_FTF_5[.,4] w_App_rate_FTF_5[.,5] ) (w_App_rate_FTF_5[.,6] w_App_rate_FTF_5[.,7] ))) ///
				(matrix(w_App_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p3) ci((w_App_rate_FTF_4[.,4] w_App_rate_FTF_4[.,5] ) (w_App_rate_FTF_4[.,6] w_App_rate_FTF_4[.,7] ))) ///		
				(matrix(w_App_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p3) ci((w_App_rate_FTF_3[.,4] w_App_rate_FTF_3[.,5] ) (w_App_rate_FTF_3[.,6] w_App_rate_FTF_3[.,7] ))) ///	
				(matrix(w_App_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p3) ci((w_App_rate_FTF_2[.,4] w_App_rate_FTF_2[.,5] ) (w_App_rate_FTF_2[.,6] w_App_rate_FTF_2[.,7] ))) ///
					, generate replace
				
				g __cgsdpos = .
				g __cgsdneg = .
				gen index = cond(!missing(__plot), _n, .)
				
				loc i=1
				display 
				foreach fert in NPS DAP App {
					foreach d in 5 4 3 2 {
						replace __cgsdpos = `w_`fert'_rate_max_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
						replace __cgsdneg = `w_`fert'_rate_min_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
						replace __cgsdpos = `w_`fert'_rate_max_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
						replace __cgsdneg = `w_`fert'_rate_min_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
					loc i = `i'+1
					}
					}

				coefplot  ///		
				(matrix(w_NPS_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((w_NPS_rate_FTF_5[.,4] w_NPS_rate_FTF_5[.,5] ) (w_NPS_rate_FTF_5[.,6] w_NPS_rate_FTF_5[.,7] ))) ///
				(matrix(w_NPS_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((w_NPS_rate_FTF_4[.,4] w_NPS_rate_FTF_4[.,5] ) (w_NPS_rate_FTF_4[.,6] w_NPS_rate_FTF_4[.,7] ))) ///		
				(matrix(w_NPS_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((w_NPS_rate_FTF_3[.,4] w_NPS_rate_FTF_3[.,5] ) (w_NPS_rate_FTF_3[.,6] w_NPS_rate_FTF_3[.,7] ))) ///	
				(matrix(w_NPS_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p1) ci((w_NPS_rate_FTF_2[.,4] w_NPS_rate_FTF_2[.,5] ) (w_NPS_rate_FTF_2[.,6] w_NPS_rate_FTF_2[.,7] ))) ///
				(matrix(w_DAP_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p2) ci((w_DAP_rate_FTF_5[.,4] w_DAP_rate_FTF_5[.,5] ) (w_DAP_rate_FTF_5[.,6] w_DAP_rate_FTF_5[.,7] ))) ///
				(matrix(w_DAP_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p2) ci((w_DAP_rate_FTF_4[.,4] w_DAP_rate_FTF_4[.,5] ) (w_DAP_rate_FTF_4[.,6] w_DAP_rate_FTF_4[.,7] ))) ///		
				(matrix(w_DAP_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p2) ci((w_DAP_rate_FTF_3[.,4] w_DAP_rate_FTF_3[.,5] ) (w_DAP_rate_FTF_3[.,6] w_DAP_rate_FTF_3[.,7] ))) ///	
				(matrix(w_DAP_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p2) ci((w_DAP_rate_FTF_2[.,4] w_DAP_rate_FTF_2[.,5] ) (w_DAP_rate_FTF_2[.,6] w_DAP_rate_FTF_2[.,7] ))) ///	
				(matrix(w_App_rate_FTF_5[,1]), label("") msymbol(O) pstyle(p3) ci((w_App_rate_FTF_5[.,4] w_App_rate_FTF_5[.,5] ) (w_App_rate_FTF_5[.,6] w_App_rate_FTF_5[.,7] ))) ///
				(matrix(w_App_rate_FTF_4[,1]), label("") msymbol(S) pstyle(p3) ci((w_App_rate_FTF_4[.,4] w_App_rate_FTF_4[.,5] ) (w_App_rate_FTF_4[.,6] w_App_rate_FTF_4[.,7] ))) ///		
				(matrix(w_App_rate_FTF_3[,1]), label("") msymbol(T) pstyle(p3) ci((w_App_rate_FTF_3[.,4] w_App_rate_FTF_3[.,5] ) (w_App_rate_FTF_3[.,6] w_App_rate_FTF_3[.,7] ))) ///	
				(matrix(w_App_rate_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p3) ci((w_App_rate_FTF_2[.,4] w_App_rate_FTF_2[.,5] ) (w_App_rate_FTF_2[.,6] w_App_rate_FTF_2[.,7] ))) ///
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(rows(6) pos(6) size(small) region(lstyle(none)) order(2 "NPS" 14 "DAP" 26 "Overall Fertilizer"  3 15 27 ///
				"5 hours" 6 18 30 "4 hours" 9 21 33 "3 hours" 12 24 36 "2 hours"  37 "" 37 "" 37 "Control group st dev"))    ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) legend(pos(6))  ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-75 75)) ylabel(-75(25)75,nogrid) ///
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))
				drop _* index
				
				graph export "Figures/Figure_S13.png", replace
			
				
	//Figure S.14: Effect of fertilizer blending facilities on yield (kg/ha): summary of estimates across different cut off points for treatment by cohort and combined across cohorts.  	
				*This code preps the control group base period st dev to be added to the coeffplot
				coefplot  ///		
				(matrix(CL_w_Yield_kgha_1_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_5[.,4] CL_w_Yield_kgha_1_FTF_5[.,5] ) (CL_w_Yield_kgha_1_FTF_5[.,6] CL_w_Yield_kgha_1_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_1_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_4[.,4] CL_w_Yield_kgha_1_FTF_4[.,5] ) (CL_w_Yield_kgha_1_FTF_4[.,6] CL_w_Yield_kgha_1_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_1_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_3[.,4] CL_w_Yield_kgha_1_FTF_3[.,5] ) (CL_w_Yield_kgha_1_FTF_3[.,6] CL_w_Yield_kgha_1_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_1_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_2[.,4] CL_w_Yield_kgha_1_FTF_2[.,5] ) (CL_w_Yield_kgha_1_FTF_2[.,6] CL_w_Yield_kgha_1_FTF_2[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_2_FTF_5[,1]), label("") msymbol(O) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_5[.,4] CL_w_Yield_kgha_2_FTF_5[.,5] ) (CL_w_Yield_kgha_2_FTF_5[.,6] CL_w_Yield_kgha_2_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_2_FTF_4[,1]), label("") msymbol(S) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_4[.,4] CL_w_Yield_kgha_2_FTF_4[.,5] ) (CL_w_Yield_kgha_2_FTF_4[.,6] CL_w_Yield_kgha_2_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_2_FTF_3[,1]), label("") msymbol(T) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_3[.,4] CL_w_Yield_kgha_2_FTF_3[.,5] ) (CL_w_Yield_kgha_2_FTF_3[.,6] CL_w_Yield_kgha_2_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_2_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_2[.,4] CL_w_Yield_kgha_2_FTF_2[.,5] ) (CL_w_Yield_kgha_2_FTF_2[.,6] CL_w_Yield_kgha_2_FTF_2[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_3_FTF_5[,1]), label("") msymbol(O) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_5[.,4] CL_w_Yield_kgha_3_FTF_5[.,5] ) (CL_w_Yield_kgha_3_FTF_5[.,6] CL_w_Yield_kgha_3_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_3_FTF_4[,1]), label("") msymbol(S) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_4[.,4] CL_w_Yield_kgha_3_FTF_4[.,5] ) (CL_w_Yield_kgha_3_FTF_4[.,6] CL_w_Yield_kgha_3_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_3_FTF_3[,1]), label("") msymbol(T) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_3[.,4] CL_w_Yield_kgha_3_FTF_3[.,5] ) (CL_w_Yield_kgha_3_FTF_3[.,6] CL_w_Yield_kgha_3_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_3_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_2[.,4] CL_w_Yield_kgha_3_FTF_2[.,5] ) (CL_w_Yield_kgha_3_FTF_2[.,6] CL_w_Yield_kgha_3_FTF_2[.,7] ))) ///
					, generate replace
				
				g __cgsdpos = .
				g __cgsdneg = .
				gen index = cond(!missing(__plot), _n, .)
				
				loc i=1
				display 
				foreach crop in 1 2 3 {
					foreach d in 5 4 3 2 {
						replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
						replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_FTF_`d'_c1' if __plot==`i' & inrange(index, 1, 24)
						replace __cgsdpos = `CL_w_Yield_kgha_`crop'_max_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
						replace __cgsdneg = `CL_w_Yield_kgha_`crop'_min_FTF_`d'_c2' if __plot==`i' & inrange(index, 25, 36)
					loc i = `i'+1
					}
					}

				coefplot  ///		
				(matrix(CL_w_Yield_kgha_1_FTF_5[,1]), label("") msymbol(O) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_5[.,4] CL_w_Yield_kgha_1_FTF_5[.,5] ) (CL_w_Yield_kgha_1_FTF_5[.,6] CL_w_Yield_kgha_1_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_1_FTF_4[,1]), label("") msymbol(S) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_4[.,4] CL_w_Yield_kgha_1_FTF_4[.,5] ) (CL_w_Yield_kgha_1_FTF_4[.,6] CL_w_Yield_kgha_1_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_1_FTF_3[,1]), label("") msymbol(T) pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_3[.,4] CL_w_Yield_kgha_1_FTF_3[.,5] ) (CL_w_Yield_kgha_1_FTF_3[.,6] CL_w_Yield_kgha_1_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_1_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p1) ci((CL_w_Yield_kgha_1_FTF_2[.,4] CL_w_Yield_kgha_1_FTF_2[.,5] ) (CL_w_Yield_kgha_1_FTF_2[.,6] CL_w_Yield_kgha_1_FTF_2[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_2_FTF_5[,1]), label("") msymbol(O) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_5[.,4] CL_w_Yield_kgha_2_FTF_5[.,5] ) (CL_w_Yield_kgha_2_FTF_5[.,6] CL_w_Yield_kgha_2_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_2_FTF_4[,1]), label("") msymbol(S) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_4[.,4] CL_w_Yield_kgha_2_FTF_4[.,5] ) (CL_w_Yield_kgha_2_FTF_4[.,6] CL_w_Yield_kgha_2_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_2_FTF_3[,1]), label("") msymbol(T) pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_3[.,4] CL_w_Yield_kgha_2_FTF_3[.,5] ) (CL_w_Yield_kgha_2_FTF_3[.,6] CL_w_Yield_kgha_2_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_2_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p2) ci((CL_w_Yield_kgha_2_FTF_2[.,4] CL_w_Yield_kgha_2_FTF_2[.,5] ) (CL_w_Yield_kgha_2_FTF_2[.,6] CL_w_Yield_kgha_2_FTF_2[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_3_FTF_5[,1]), label("") msymbol(O) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_5[.,4] CL_w_Yield_kgha_3_FTF_5[.,5] ) (CL_w_Yield_kgha_3_FTF_5[.,6] CL_w_Yield_kgha_3_FTF_5[.,7] ))) ///
				(matrix(CL_w_Yield_kgha_3_FTF_4[,1]), label("") msymbol(S) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_4[.,4] CL_w_Yield_kgha_3_FTF_4[.,5] ) (CL_w_Yield_kgha_3_FTF_4[.,6] CL_w_Yield_kgha_3_FTF_4[.,7] ))) ///		
				(matrix(CL_w_Yield_kgha_3_FTF_3[,1]), label("") msymbol(T) pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_3[.,4] CL_w_Yield_kgha_3_FTF_3[.,5] ) (CL_w_Yield_kgha_3_FTF_3[.,6] CL_w_Yield_kgha_3_FTF_3[.,7] ))) ///	
				(matrix(CL_w_Yield_kgha_3_FTF_2[,1]), label("") msymbol(X) msize(large)  pstyle(p3) ci((CL_w_Yield_kgha_3_FTF_2[.,4] CL_w_Yield_kgha_3_FTF_2[.,5] ) (CL_w_Yield_kgha_3_FTF_2[.,6] CL_w_Yield_kgha_3_FTF_2[.,7] ))) ///	
				, rename (r1= Combined r2 = "Cohort one" r3="Cohort two")  ///
				drop (Post Treat *.Treat *.Post  *.Woreda_U_id_code land_area_sqm_10thou Female Litrate _cons Pre_avg Post_avg) ///
				vertical legend(rows(6) pos(6) size(small) region(lstyle(none)) order(2 "Teff" 14 "Wheat" 26 "Maize"  3 15 27 ///
				"5 hours" 6 18 30 "4 hours" 9 21 33 "3 hours" 12 24 36 "2 hours"  37 "" 37 "" 37 "Control group st dev"))    ///
				graphregion(color(white)) coeflabels( Treatment="Treatment Effects Coefficients")  ///
				msymbol(D)  levels( 95 90) legend(pos(6))  ///
				yline(0, lstyle(1) lcolor(red)) yscale(r(-1200 1200)) ylabel(-1000(500)1000,nogrid)	///	
				addplot((scatter __cgsdpos __at, msymbol(smplus) mcolor(gs8) msize(medium)) (scatter __cgsdneg __at, msymbol(smplus) mcolor(gs8) msize(medium)))
				drop _* index
				
				graph export "Figures/Figure_S14.png", replace
	
	
	

* Tables - mian manuscript  
				use FTF.dta, clear 
					
	//Table 1: Average values of pre-treatment covariates and outcome variables by treatment status and cohort
				*Cohort 1 pre
				iebaltab w_area_planted Female Litrate  NPS_Dummy DAP_Dummy Fert_Dummy w_NPS_rate w_DAP_rate w_App_rate ///
				if year==1, grpvar(Treatment_1) stats(desc(sd)) ///
				savexlsx ("Tables\Summary_Tr_Cn_c1_pre.xlsx") rowvarlabels rowlabels() ///
				order (0 1) grpl(0 Control @ 1 Treated)  tot replace 

				*Cohort 2 pre
				iebaltab w_area_planted Female Litrate NPS_Dummy DAP_Dummy Fert_Dummy w_NPS_rate w_DAP_rate w_App_rate ///
				if year==2 & first_treat!=2, grpvar(Treatment_2) stats(desc(sd)) ///
				savexlsx ("Tables\Summary_Tr_Cn_c2_pre.xlsx") rowvarlabels rowlabels() ///
				order (0 1) grpl(0 Control @ 1 Treated)  tot replace  
				
				*Cohort 1 post
				iebaltab NPS_Dummy DAP_Dummy Fert_Dummy w_NPS_rate w_DAP_rate w_App_rate ///
				if year==2, grpvar(Treatment_1) stats(desc(sd)) ///
				savexlsx ("Tables\Summary_Tr_Cn_c1_post.xlsx") rowvarlabels rowlabels() ///
				order (0 1) grpl(0 Control @ 1 Treated)  tot replace 
				
				*Cohort 2 post
				iebaltab NPS_Dummy DAP_Dummy Fert_Dummy w_NPS_rate w_DAP_rate w_App_rate ///
				if year==3 & first_treat!=2, grpvar(Treatment_2) stats(desc(sd)) ///
				savexlsx ("Tables\Summary_Tr_Cn_c2_post.xlsx") rowvarlabels rowlabels() ///
				order (0 1) grpl(0 Control @ 1 Treated)  tot replace 
		
	//Table 2: Correspondence between datasets and study time periods. 
				*-This table is manually created without a code that can replicate it.
			
	//Table 3: Average fertilizer prices (birr/kg) by type, over time, across groups of farmers 
				gen NPS_DAP = NPS_Price_mr/DAP_Price_mr
				gen DAP1_tempo =  DAP_Price_mr if year ==1
				gen DAP3_tempo =  DAP_Price_mr if year ==3
				bysort id: egen DAP1 = max(DAP1_tempo)
				bysort id: egen DAP3 = max(DAP3_tempo)
				gen DAP3_DAP1 = DAP3/DAP1 if year == 3
			
				*Never treated 
				estpost tabstat NPS_Price_mr DAP_Price_mr NPS_DAP DAP3_DAP1   if first_treat==0 /*Third cohort*/, by(year ) statistics( mean sd ) columns(statistics) nototal
				esttab using "Tables\Table_3.rtf", main(mean  %9.2f sd %9.2f min %9.2f max %9.2f ) aux(sd  %9.2f min %9.2f  max %9.2f )    ///
				title("Table##: Average fertilizer prices (birr/kg) by type, over time, across groups of farmers -Never Treated  ") unstack   nonote obslast nostar nogaps label rtf replace  
				*First cohort
				estpost tabstat NPS_Price_mr DAP_Price_mr NPS_DAP DAP3_DAP1   if first_treat==2 /*First cohort*/, by(year ) statistics( mean sd ) columns(statistics) nototal
				esttab using "Tables\Table_3.rtf", main(mean  %9.2f sd %9.2f min %9.2f max %9.2f ) aux(sd  %9.2f min %9.2f  max %9.2f )    ///
				title("Table##: Average fertilizer prices (birr/kg) by type, over time, across groups of farmers - First Cohort ") unstack nonote obslast nostar nogaps label rtf append  
				*Second Cohort
				estpost tabstat NPS_Price_mr DAP_Price_mr NPS_DAP DAP3_DAP1   if first_treat==3 /*First cohort*/, by(year ) statistics( mean sd ) columns(statistics) nototal
				esttab using "Tables\Table_3.rtf", main(mean  %9.2f sd %9.2f min %9.2f max %9.2f ) aux(sd  %9.2f min %9.2f  max %9.2f )    ///
				title("Table##: Average fertilizer prices (birr/kg) by type, over time, across groups of farmers  - Second treated  ") unstack nonote obslast nostar nogaps label rtf append  
		
	//Table 4: Cross tabulation of household level observations by fertilizer blending facility and demonstration plot treatment status.
				*Table 4 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.
		
	//Table 5: Heterogeneity in the impact of FBFS on fertilizer adoption, application rates, and yields over access to demonstration plots.
				*Table 5 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.	
		
	//Table 6: Pre-trend hypothesis tests for core outcomes. 
				*Table 6 is constructed using the Agricultural Sample Survey (AgSS) data. Due to data access restrictions imposed by the data provider, we are unable to make the data publicly available.

				
				
* Tables - Supplemental Materials  
		
	//Table S.1: Number of DPs with NPS demonstration plots by region, crop, and type of demonstration plot during the 2013/14 cropping year.
				*-This table is manually created without a code that can replicate it.
			
	//Table S.2 - S.19 
				*this is just working on the table code -- to be pasted into the main do file.
				** This version -- use the csdid actual specification


				adopath + "ado/"

				loc seed 19810310
				loc covars land_area_sqm_10thou Female Litrate 
				loc ncovar: word count `covars'

				loc clustsizelist A05 A04 A06
				loc distlist 100 200 300
				loc siglist 0.1 0.05 0.01

				loc clustsize_default A05
				loc clustsize_alt A04
				loc clustsize_spatial A05 /*A25 (ea) and A04 (woreda) don't work. All units in the spatial clusters are assumed to have the same location.*/

				use FTF.dta, clear 
				cap g A04 = Woreda_U_id_code
				cap g lat = automatic_ /*This is in decimal degrees*/
				cap g lon = automati_1 /*This is in decimal degrees*/
				bysort `clustsize_spatial': egen s_1 = median(lat)
				bysort `clustsize_spatial': egen s_2 = median(lon) 
				*scatter s_1 s_2

				* Rename outcome variables so they are short enough for code

				loc NPS_Dummy_label "NPS adoption (binary)"
				loc NPS_Dummy_Label "NPS adoption (binary)"
				loc DAP_Dummy_label "DAP adoption (binary)"
				loc DAP_Dummy_Label "DAP adoption (binary)"
				rename Fert_Dummy F_Dummy /*variable name needs to be short for the code to run*/
				loc F_Dummy_label "any fertilizer adoption (binary)"
				loc F_Dummy_Label "Any fert adoption (binary)"
				rename w_NPS_rate NPS_Cont 
				loc NPS_Cont_label "NPS application rate (kg/ha)"
				loc NPS_Cont_Label "NPS application rate (kg/ha)"
				rename w_DAP_rate DAP_Cont 
				loc DAP_Cont_label "DAP application rate (kg/ha)"
				loc DAP_Cont_Label "DAP application rate (kg/ha)"
				rename w_App_rate F_Cont 
				loc F_Cont_label "fertilizer application rate (kg/ha)"
				loc F_Cont_Label "Fert application rate (kg/ha)"
				rename CL_w_Yield_kgha_1 Yield_1 
				loc Yield_1_label "teff yields (kg/ha)"
				loc Yield_1_Label "Teff yields (kg/ha)"
				rename CL_w_Yield_kgha_2 Yield_2 
				loc Yield_2_label "wheat yields (kg/ha)"
				loc Yield_2_Label "Wheat yields (kg/ha)"
				rename CL_w_Yield_kgha_3 Yield_3
				loc Yield_3_label "maize yields (kg/ha)"
				loc Yield_3_Label "Maize yields (kg/ha)"
				rename CL_w_area_1 Area_1 
				loc Area_1_label "teff area planted (ha)"
				loc Area_1_Label "Teff area planted (ha)"
				rename CL_w_area_2 Area_2 
				loc Area_2_label "wheat area planted (ha)"
				loc Area_2_Label "Wheat area planted (ha)"
				rename CL_w_area_3 Area_3 
				loc Area_3_label "maize area planted (ha)"
				loc Area_3_Label "Maize area planted (ha)"
				rename w_area_planted Area_t
				loc Area_t_label "total area planted (ha)"
				loc Area_t_Label "Total area planted (ha)"
				rename TOT_V_harvest_kg v_harv
				loc v_harv_label "total harvest value (ETB)"
				loc v_harv_Label "Total harvest value (ETB)"
				rename CL_V_w_harvest_kg_1 v_harv1
				loc v_harv1_label "total teff harvest value (ETB)"
				loc v_harv1_Label "Total teff harvest value (ETB)"
				rename CL_V_w_harvest_kg_2 v_harv2
				loc v_harv2_label "total wheat harvest value (ETB)"
				loc v_harv2_Label "Total wheat harvest value (ETB)"
				rename CL_V_w_harvest_kg_3 v_harv3
				loc v_harv3_label "total maize harvest value (ETB)"
				loc v_harv3_Label "Total maize harvest value (ETB)"
				rename DAP_Price_mr DAP_pr 
				loc DAP_pr_label "DAP price (birr/kg)"
				loc DAP_pr_Label "DAP price (birr/kg)"
				rename pcexp_aeu pcexpae
				loc pcexpae_label "expenditures per adult equiv (ETB/day)" /*CONFIRM*/
				loc pcexpae_Label "Expenditures per adult equiv (ETB/day)" /*CONFIRM*/
				rename paexpD_R pcexpaeDR
				loc pcexpaeDR_label "real expenditures per adult equivalent (2011 ETB/day)" 
				loc pcexpaeDR_Label "Exp per adult equiv (2011 ETB/day)" 

				la var land_area_sqm_10thou "Landholdings (ha)"
				la var Female "Female headed"
				la var Litrate "Literate (head)"


				forvalues c = 1/3 {
				g Y`c'm = Yield_`c'==. 
				bysort id: egen Y`c'm_ever = max(Y`c'm)
				replace Yield_`c' = . if Y`c'm_ever==1

				g vh`c'm = v_harv`c'==.
				bysort id: egen vh`c'm_ever = max(vh`c'm)
				replace v_harv`c' = . if vh`c'm_ever==1	
				}

				loc outcomelist NPS_Dummy DAP_Dummy F_Dummy NPS_Cont DAP_Cont F_Cont DAP_pr Yield_1 Yield_2 Yield_3 v_harv v_harv1 v_harv2 v_harv3  pcexpaeDR


				xtset


				** First Cohort periods 1 & 2 DID setup
				g cohort1atreat = (first_treat==2)
				g cohort1acontrol = (first_treat==3 | first_treat==0)
				g cohort1apre = (year==1)
				g cohort1apost = (year==2)
				g cohort1atreatpost = cohort1atreat * cohort1apost
				g cohort1atoinclude = (cohort1atreat==1 | cohort1acontrol==1) & (cohort1apre==1 | cohort1apost==1)
				bysort cohort1atoinclude: tab first_treat year

				** First Cohort periods 1 & 3 DID setup
				g cohort1btreat = (first_treat==2)
				g cohort1bcontrol = (first_treat==0)
				g cohort1bpre = (year==1)
				g cohort1bpost = (year==3)
				g cohort1btreatpost = cohort1btreat * cohort1bpost
				g cohort1btoinclude = (cohort1btreat==1 | cohort1bcontrol==1) & (cohort1bpre==1 | cohort1bpost==1)
				bysort cohort1btoinclude: tab first_treat year

				** Second period DID setup
				g cohort2treat = (first_treat==3)
				g cohort2control = (first_treat==0)
				g cohort2pre = (year==2)
				g cohort2post = (year==3)
				g cohort2treatpost = cohort2treat * cohort2post
				g cohort2toinclude = (cohort2treat==1 | cohort2control==1) & (cohort2pre==1 | cohort2post==1)
				bysort cohort2toinclude: tab first_treat year

				foreach c in 1a 1b 2 {
				la var cohort`c'treat "Treat"
				la var cohort`c'post "Post"
				la var cohort`c'treatpost "Treat X Post"
				}

				** Create a covariate vector for periods 1a and 1b period 2
				forvalues i = 1/2 {
				loc covars`i' ""
				foreach var in `covars'{
					tempvar var`i'temp 
					g `var`i'temp' = `var' if year==`i'
					bysort id: egen `var'`i' = median(`var`i'temp')
					loc covars`i' `covars`i'' `var'`i'
					local lbl : variable label `var' 
					label var `var'`i' "`lbl'"
					}
				} 

				display "`covars1'"
				display "`covars2'"
				loc covars1b `covars1'
				loc covars1a `covars1'

				loc interactlist 
				g y1=(year==1)
				g y2=(year==2)
				g y3=(year==3)
				g ft2=(first_treat==2)
				g ft3=(first_treat==3)
				g ft0=(first_treat==0)
				g y1ft2 = y1*ft2
				g y1ft3 = y1*ft3
				g y1ft0 = y1*ft0
				g y2ft2 = y2*ft2
				g y2ft3 = y2*ft3
				g y2ft0 = y2*ft0
				g y3ft2 = y3*ft2
				g y3ft3 = y3*ft3
				g y3ft0 = y3*ft0
				loc ylist y1 y2 y3
				loc ftlist ft2 ft3 ft0
				loc yftlist y1ft2 y1ft3 y1ft0 y2ft2 y2ft3 y2ft0 y3ft2 y3ft3 y3ft0 

				loc i = 1
				loc yftclist  ""
				foreach covar in `covars1' {
				foreach var in `yftlist' {
					g `var'c`i' = `var'*`covar'
					loc yftclist `yftclist' `var'c`i'
					}
				loc i = `i'+1
				}
				loc xvarlist `ylist' `ftlist' `yftlist' `yftclist' /*`covars1'*/


				** Save analysis dataset

				tempfile dataset_all
				save `dataset_all', replace

				foreach c in 1a 1b 2 {
				use `dataset_all', clear
				keep if cohort`c'toinclude==1
				tempfile dataset_c`c'
				save `dataset_c`c'', replace
				}


				** CSDID results for main tables
				*
				foreach outcome in `outcomelist' {
				use `dataset_all', clear
				loc coeflist ""

				* Get cluster size and N by cohort and use naming for the stars table
				foreach clust in `clustsizelist' {
					unique `clust' if (`outcome'!=.) & (cohort1atoinclude==1 | cohort1btoinclude==1 | cohort2toinclude==1)
					loc `outcome'_coef_csdidga_NC`clust' = r(unique)
					loc `outcome'_coef_csdidga_N`clust' = r(N)
					
					unique `clust' if (`outcome'!=.) & (cohort1atoinclude==1 | cohort1btoinclude==1)
					loc `outcome'_coef_csdidg2_NC`clust' = r(unique)
					loc `outcome'_coef_csdidg2_N`clust' = r(N)

					unique `clust' if (`outcome'!=.) & (cohort2toinclude==1)
					loc `outcome'_coef_csdidg3_NC`clust' = r(unique)
					loc `outcome'_coef_csdidg3_N`clust' = r(N)

					unique `clust' if (`outcome'!=.) & (cohort1atoinclude==1)
					loc `outcome'_coef_csdidc1a_NC`clust' = r(unique)
					loc `outcome'_coef_csdidc1a_N`clust' = r(N)

					unique `clust' if (`outcome'!=.) & (cohort1btoinclude==1)
					loc `outcome'_coef_csdidc1b_NC`clust' = r(unique)
					loc `outcome'_coef_csdidc1b_N`clust' = r(N)
					
					unique `clust' if (`outcome'!=.) & (cohort2toinclude==1)
					loc `outcome'_coef_csdidc2_NC`clust' = r(unique)
					loc `outcome'_coef_csdidc2_N`clust' = r(N)
					}

					
				*By group
				csdid `outcome' `covars', ivar(id) time(year) gvar(first_treat) method(dripw) notyet agg(group) wboot cluster(`clustsize_default') rseed(`seed') level(95) 
				estimates store csdid_`outcome'_byg
				loc `outcome'_N_csdid = e(N)
				mat coef=e(b)
				loc `outcome'_coef_csdidga = coef[1,1]
				loc `outcome'_coef_csdidg2 = coef[1,2]
				loc `outcome'_coef_csdidg3 = coef[1,3]
					
				su `outcome' if (cohort1atoinclude==1 | cohort1btoinclude==1) & year==1 & (cohort1acontrol==1 | cohort1bcontrol==1)
				loc `outcome'_cmean_ga = r(mean)
				loc `outcome'_csd_ga = r(sd)
				loc `outcome'_cmean_g2 = r(mean)
				loc `outcome'_csd_g2 = r(sd)
				su `outcome' if (cohort2toinclude==1) & year==2 & (cohort2control==1)
				loc `outcome'_cmean_g3 = r(mean) 
				loc `outcome'_csd_g3 = r(sd) 

				foreach g in a 2 3 {
					loc `outcome'_coef_csdidg`g' : display %09.3fc ``outcome'_coef_csdidg`g''
					loc `outcome'_coef_csdidg`g'_N : display %09.0fc ``outcome'_N_csdid'
					loc coeflist `coeflist' `outcome'_coef_csdidg`g'
					loc `outcome'_cmean_g`g' : display %09.3fc ``outcome'_cmean_g`g''
					loc `outcome'_csd_g`g' : display %09.3fc ``outcome'_csd_g`g''
					}
					
				*By gt
				csdid `outcome' `covars', ivar(id) time(year) gvar(first_treat) method(dripw) notyet agg(attgt) wboot cluster(`clustsize_default') rseed(`seed') level(95) 
				estimates store csdid_`outcome'_bygt
				mat coef=e(b)
				mat N=e(gtt)
				loc `outcome'_coef_csdidc1a = coef[1,1]
				loc `outcome'_N_csdidc1a = N[1,4]
				loc `outcome'_coef_csdidc1b = coef[1,2]
				loc `outcome'_N_csdidc1b = N[2,4]
				loc `outcome'_coef_csdidc2 = coef[1,4]
				loc `outcome'_N_csdidc2 = N[4,4]

				su `outcome' if (cohort1atoinclude==1) & year==1 & (cohort1acontrol==1)
				loc `outcome'_cmean_c1a = r(mean)
				loc `outcome'_csd_c1a = r(sd)
				su `outcome' if (cohort1btoinclude==1) & year==1 & (cohort1bcontrol==1)
				loc `outcome'_cmean_c1b = r(mean)
				loc `outcome'_csd_c1b = r(sd)
				su `outcome' if (cohort2toinclude==1) & year==2 & (cohort2control==1)
				loc `outcome'_cmean_c2 = r(mean) 
				loc `outcome'_csd_c2 = r(sd) 

				foreach c in 1a 1b 2 {
					loc `outcome'_coef_csdidc`c' : display %09.3fc ``outcome'_coef_csdidc`c''
					loc coeflist `coeflist' `outcome'_coef_csdidc`c'
					loc `outcome'_N_csdidc`c' : display %09.0fc ``outcome'_N_csdidc`c''
					loc `outcome'_cmean_c`c' : display %09.3fc ``outcome'_cmean_c`c''
					loc `outcome'_csd_c`c' : display %09.3fc ``outcome'_csd_c`c''
					}

				foreach clust in `clustsizelist' {
					*By group
					csdid `outcome' `covars', ivar(id) time(year) gvar(first_treat) method(dripw) notyet agg(group) wboot cluster(`clust') rseed(`seed') level(95) 
					estimates store csdid_`outcome'_byg_`clust'
					mat var = e(V)
					loc `outcome'_coef_csdidga_se_`clust' = sqrt(var[1,1])
					loc `outcome'_coef_csdidg2_se_`clust' = sqrt(var[2,2])
					loc `outcome'_coef_csdidg3_se_`clust' = sqrt(var[3,3])
					
					loc `outcome'_Nclust_se_`clust' = e(N_clust)
					loc `outcome'_Nclust_se_`clust' : display %09.0fc ``outcome'_Nclust_se_`clust''
					
					foreach g in a 2 3 {
						loc `outcome'_coef_csdidg`g'_se_`clust' : display %09.3fc ``outcome'_coef_csdidg`g'_se_`clust''
						loc `outcome'_Nclust_g`g'_se_`clust' : display %09.0fc ``outcome'_Nclust_se_`clust''
						}
					
					*By gt
					csdid `outcome' `covars', ivar(id) time(year) gvar(first_treat) method(dripw) notyet agg(attgt) wboot cluster(`clust') rseed(`seed') level(95) 
					estimates store csdid_`outcome'_bygt_`clust'
					mat var = e(V)
					loc `outcome'_coef_csdidc1a_se_`clust' = sqrt(var[1,1])
					loc `outcome'_coef_csdidc1b_se_`clust' = sqrt(var[2,2])
					loc `outcome'_coef_csdidc2_se_`clust' = sqrt(var[4,4])
					foreach c in 1a 1b 2 {
						loc `outcome'_coef_csdidc`c'_se_`clust' : display %09.3fc ``outcome'_coef_csdidc`c'_se_`clust''
						}
					
					display "`coeflist'"
					foreach coef in `coeflist' {
						loc `coef'_st_`clust' = ""
						display "`coef'_st_`clust'"
						display "`coef'_se_`clust'"
						foreach sig in `siglist' {
							loc tstat = 0 /*default tstat*/
							*cap loc tstat = (abs(``coef''/``coef'_se_`clust'')>invttail(``outcome'_Nclust_se_`clust'',`sig'/2))
							cap loc tstat = (abs(``coef''/``coef'_se_`clust'')>invttail(``coef'_NC`clust'',`sig'/2))
							if `tstat'==float(1) {
								loc `coef'_st_`clust' ``coef'_st_`clust''*  
								}
							display "coef: `coef'; sig: `sig' tstat: `tstat'"
							}
						display "coef: `coef'; stars: ``coef'_st_`clust''"	
						}
					}		
				}


				** Make the 3 main results tables: Table S.2, S.3 and S.4

				loc varlistcsdidtable1 NPS_Dummy DAP_Dummy F_Dummy NPS_Cont DAP_Cont F_Cont
				loc varlistcsdidtable2 Yield_1 Yield_2 Yield_3 /*Area_t*/ DAP_pr
				loc varlistcsdidtable3 v_harv v_harv1 v_harv2 v_harv3 pcexpaeDR

	//Table S.2: Effect of FBFs on Fertilizer Adoption Overall and by cohort - Doubly Robust Difference in Difference Model.
				cap texdoc close
				texdoc init "Tables/FertDummyCont_main_csdid", replace force 
				tex \begin{table}
				tex \begin{adjustbox}{width=.8\linewidth, keepaspectratio}
				tex \begin{threeparttable}
				tex \caption{Effect of FBFs on Fertilizer Adoption Overall and by cohort - Doubly Robust Difference in Difference Model}\label{results_csdid_FertDummyCont}
				tex \begin{tabular}{lclclclclclcl}
				tex \hline \hline
				tex & \multicolumn{2}{c}{(1)} & \multicolumn{2}{c}{(2)} & \multicolumn{2}{c}{(3)} & \multicolumn{2}{c}{(4)} & \multicolumn{2}{c}{(5)} & \multicolumn{2}{c}{(6)}\\
				tex & \multicolumn{2}{c}{NPS} & \multicolumn{2}{c}{DAP} & \multicolumn{2}{c}{Fert} & \multicolumn{2}{c}{NPS} & \multicolumn{2}{c}{DAP} & \multicolumn{2}{c}{Fert}\\
				tex & \multicolumn{2}{c}{(binary)} & \multicolumn{2}{c}{(binary)} & \multicolumn{2}{c}{(binary)} & \multicolumn{2}{c}{(kg/ha)} & \multicolumn{2}{c}{(kg/ha)} & \multicolumn{2}{c}{(kg/ha)}\\
				tex \hline
				tex Overall & `NPS_Dummy_coef_csdidga' &  & `DAP_Dummy_coef_csdidga' & & `F_Dummy_coef_csdidga' & & `NPS_Cont_coef_csdidga' & & `DAP_Cont_coef_csdidga' & & `F_Cont_coef_csdidga' & \\
				tex & (`NPS_Dummy_coef_csdidga_se_A05') &  `NPS_Dummy_coef_csdidga_st_A05' & (`DAP_Dummy_coef_csdidga_se_A05') &  `DAP_Dummy_coef_csdidga_st_A05' & (`F_Dummy_coef_csdidga_se_A05') &  `F_Dummy_coef_csdidga_st_A05' & (`NPS_Cont_coef_csdidga_se_A05') &  `NPS_Cont_coef_csdidga_st_A05' & (`DAP_Cont_coef_csdidga_se_A05') &  `DAP_Cont_coef_csdidga_st_A05' & (`F_Cont_coef_csdidga_se_A05') &  `F_Cont_coef_csdidga_st_A05' \\
				tex & [`NPS_Dummy_coef_csdidga_se_A04'] &  `NPS_Dummy_coef_csdidga_st_A04' & [`DAP_Dummy_coef_csdidga_se_A04'] &  `DAP_Dummy_coef_csdidga_st_A04' & [`F_Dummy_coef_csdidga_se_A04'] &  `F_Dummy_coef_csdidga_st_A04' & [`NPS_Cont_coef_csdidga_se_A04'] &  `NPS_Cont_coef_csdidga_st_A04' & [`DAP_Cont_coef_csdidga_se_A04'] &  `DAP_Cont_coef_csdidga_st_A04' & [`F_Cont_coef_csdidga_se_A04'] &  `F_Cont_coef_csdidga_st_A04' \\
				tex Cohort 1a & `NPS_Dummy_coef_csdidc1a' &  & `DAP_Dummy_coef_csdidc1a' & & `F_Dummy_coef_csdidc1a' & & `NPS_Cont_coef_csdidc1a' & & `DAP_Cont_coef_csdidc1a' & & `F_Cont_coef_csdidc1a' & \\
				tex & (`NPS_Dummy_coef_csdidc1a_se_A05') &  `NPS_Dummy_coef_csdidc1a_st_A05' & (`DAP_Dummy_coef_csdidc1a_se_A05') &  `DAP_Dummy_coef_csdidc1a_st_A05' & (`F_Dummy_coef_csdidc1a_se_A05') &  `F_Dummy_coef_csdidc1a_st_A05' & (`NPS_Cont_coef_csdidc1a_se_A05') &  `NPS_Cont_coef_csdidc1a_st_A05' & (`DAP_Cont_coef_csdidc1a_se_A05') &  `DAP_Cont_coef_csdidc1a_st_A05' & (`F_Cont_coef_csdidc1a_se_A05') &  `F_Cont_coef_csdidc1a_st_A05' \\
				tex & [`NPS_Dummy_coef_csdidc1a_se_A04'] &  `NPS_Dummy_coef_csdidc1a_st_A04' & [`DAP_Dummy_coef_csdidc1a_se_A04'] &  `DAP_Dummy_coef_csdidc1a_st_A04' & [`F_Dummy_coef_csdidc1a_se_A04'] &  `F_Dummy_coef_csdidc1a_st_A04' & [`NPS_Cont_coef_csdidc1a_se_A04'] &  `NPS_Cont_coef_csdidc1a_st_A04' & [`DAP_Cont_coef_csdidc1a_se_A04'] &  `DAP_Cont_coef_csdidc1a_st_A04' & [`F_Cont_coef_csdidc1a_se_A04'] &  `F_Cont_coef_csdidc1a_st_A04' \\
				tex Cohort 1b & `NPS_Dummy_coef_csdidc1b' &  & `DAP_Dummy_coef_csdidc1b' & & `F_Dummy_coef_csdidc1b' & & `NPS_Cont_coef_csdidc1b' & & `DAP_Cont_coef_csdidc1b' & & `F_Cont_coef_csdidc1b' & \\
				tex & (`NPS_Dummy_coef_csdidc1b_se_A05') &  `NPS_Dummy_coef_csdidc1b_st_A05' & (`DAP_Dummy_coef_csdidc1b_se_A05') &  `DAP_Dummy_coef_csdidc1b_st_A05' & (`F_Dummy_coef_csdidc1b_se_A05') &  `F_Dummy_coef_csdidc1b_st_A05' & (`NPS_Cont_coef_csdidc1b_se_A05') &  `NPS_Cont_coef_csdidc1b_st_A05' & (`DAP_Cont_coef_csdidc1b_se_A05') &  `DAP_Cont_coef_csdidc1b_st_A05' & (`F_Cont_coef_csdidc1b_se_A05') &  `F_Cont_coef_csdidc1b_st_A05' \\
				tex & [`NPS_Dummy_coef_csdidc1b_se_A04'] &  `NPS_Dummy_coef_csdidc1b_st_A04' & [`DAP_Dummy_coef_csdidc1b_se_A04'] &  `DAP_Dummy_coef_csdidc1b_st_A04' & [`F_Dummy_coef_csdidc1b_se_A04'] &  `F_Dummy_coef_csdidc1b_st_A04' & [`NPS_Cont_coef_csdidc1b_se_A04'] &  `NPS_Cont_coef_csdidc1b_st_A04' & [`DAP_Cont_coef_csdidc1b_se_A04'] &  `DAP_Cont_coef_csdidc1b_st_A04' & [`F_Cont_coef_csdidc1b_se_A04'] &  `F_Cont_coef_csdidc1b_st_A04' \\
				tex Cohort 1 combined & `NPS_Dummy_coef_csdidg2' &  & `DAP_Dummy_coef_csdidg2' & & `F_Dummy_coef_csdidg2' & & `NPS_Cont_coef_csdidg2' & & `DAP_Cont_coef_csdidg2' & & `F_Cont_coef_csdidg2' & \\
				tex & (`NPS_Dummy_coef_csdidg2_se_A05') &  `NPS_Dummy_coef_csdidg2_st_A05' & (`DAP_Dummy_coef_csdidg2_se_A05') &  `DAP_Dummy_coef_csdidg2_st_A05' & (`F_Dummy_coef_csdidg2_se_A05') &  `F_Dummy_coef_csdidg2_st_A05' & (`NPS_Cont_coef_csdidg2_se_A05') &  `NPS_Cont_coef_csdidg2_st_A05' & (`DAP_Cont_coef_csdidg2_se_A05') &  `DAP_Cont_coef_csdidg2_st_A05' & (`F_Cont_coef_csdidg2_se_A05') &  `F_Cont_coef_csdidg2_st_A05' \\
				tex & [`NPS_Dummy_coef_csdidg2_se_A04'] &  `NPS_Dummy_coef_csdidg2_st_A04' & [`DAP_Dummy_coef_csdidg2_se_A04'] &  `DAP_Dummy_coef_csdidg2_st_A04' & [`F_Dummy_coef_csdidg2_se_A04'] &  `F_Dummy_coef_csdidg2_st_A04' & [`NPS_Cont_coef_csdidg2_se_A04'] &  `NPS_Cont_coef_csdidg2_st_A04' & [`DAP_Cont_coef_csdidg2_se_A04'] &  `DAP_Cont_coef_csdidg2_st_A04' & [`F_Cont_coef_csdidg2_se_A04'] &  `F_Cont_coef_csdidg2_st_A04' \\
				tex Cohort 2 & `NPS_Dummy_coef_csdidg3' &  & `DAP_Dummy_coef_csdidg3' & & `F_Dummy_coef_csdidg3' & & `NPS_Cont_coef_csdidg3' & & `DAP_Cont_coef_csdidg3' & & `F_Cont_coef_csdidg3' & \\
				tex & (`NPS_Dummy_coef_csdidg3_se_A05') &  `NPS_Dummy_coef_csdidg3_st_A05' & (`DAP_Dummy_coef_csdidg3_se_A05') &  `DAP_Dummy_coef_csdidg3_st_A05' & (`F_Dummy_coef_csdidg3_se_A05') &  `F_Dummy_coef_csdidg3_st_A05' & (`NPS_Cont_coef_csdidg3_se_A05') &  `NPS_Cont_coef_csdidg3_st_A05' & (`DAP_Cont_coef_csdidg3_se_A05') &  `DAP_Cont_coef_csdidg3_st_A05' & (`F_Cont_coef_csdidg3_se_A05') &  `F_Cont_coef_csdidg3_st_A05' \\
				tex & [`NPS_Dummy_coef_csdidg3_se_A04'] &  `NPS_Dummy_coef_csdidg3_st_A04' & [`DAP_Dummy_coef_csdidg3_se_A04'] &  `DAP_Dummy_coef_csdidg3_st_A04' & [`F_Dummy_coef_csdidg3_se_A04'] &  `F_Dummy_coef_csdidg3_st_A04' & [`NPS_Cont_coef_csdidg3_se_A04'] &  `NPS_Cont_coef_csdidg3_st_A04' & [`DAP_Cont_coef_csdidg3_se_A04'] &  `DAP_Cont_coef_csdidg3_st_A04' & [`F_Cont_coef_csdidg3_se_A04'] &  `F_Cont_coef_csdidg3_st_A04' \\ 
				tex \hline 
				tex Pre-period control mean (sd) & & & & & & & & & & & & \\
				tex \hspace{3mm} Overall & `NPS_Dummy_cmean_ga' & & `DAP_Dummy_cmean_ga' & & `F_Dummy_cmean_ga' & & `NPS_Cont_cmean_ga' & & `DAP_Cont_cmean_ga' & & `F_Cont_cmean_ga' & \\
				tex  & (`NPS_Dummy_csd_ga') & & (`DAP_Dummy_csd_ga') & & (`F_Dummy_csd_ga') & & (`NPS_Cont_csd_ga') & & (`DAP_Cont_csd_ga') & & (`F_Cont_csd_ga') & \\
				tex \hspace{3mm} Cohort 1a & `NPS_Dummy_cmean_c1a' & & `DAP_Dummy_cmean_c1a' & & `F_Dummy_cmean_c1a' & & `NPS_Cont_cmean_c1a' & & `DAP_Cont_cmean_c1a' & & `F_Cont_cmean_c1a' & \\
				tex & (`NPS_Dummy_csd_c1a') & & (`DAP_Dummy_csd_c1a') & & (`F_Dummy_csd_c1a') & & (`NPS_Cont_csd_c1a') & & (`DAP_Cont_csd_c1a') & & (`F_Cont_csd_c1a') & \\
				tex \hspace{3mm} Cohort 1b & `NPS_Dummy_cmean_c1b' & & `DAP_Dummy_cmean_c1b' & & `F_Dummy_cmean_c1b' & & `NPS_Cont_cmean_c1b' & & `DAP_Cont_cmean_c1b' & & `F_Cont_cmean_c1b' & \\
				tex & (`NPS_Dummy_csd_c1b') & & (`DAP_Dummy_csd_c1b') & & (`F_Dummy_csd_c1b') & & (`NPS_Cont_csd_c1b') & & (`DAP_Cont_csd_c1b') & & (`F_Cont_csd_c1b') & \\
				tex \hspace{3mm} Cohort 1 combined & `NPS_Dummy_cmean_g2' & & `DAP_Dummy_cmean_g2' & & `F_Dummy_cmean_g2' & & `NPS_Cont_cmean_g2' & & `DAP_Cont_cmean_g2' & & `F_Cont_cmean_g2' & \\
				tex  & (`NPS_Dummy_csd_g2') & & (`DAP_Dummy_csd_g2') & & (`F_Dummy_csd_g2') & & (`NPS_Cont_csd_g2') & & (`DAP_Cont_csd_g2') & & (`F_Cont_csd_g2') & \\
				tex \hspace{3mm} Cohort 2 & `NPS_Dummy_cmean_g3' & & `DAP_Dummy_cmean_g3' & & `F_Dummy_cmean_g3' & & `NPS_Cont_cmean_g3' & & `DAP_Cont_cmean_g3' & & `F_Cont_cmean_g3' & \\
				tex  & (`NPS_Dummy_csd_g3') & & (`DAP_Dummy_csd_g3') & & (`F_Dummy_csd_g3') & & (`NPS_Cont_csd_g3') & & (`DAP_Cont_csd_g3') & & (`F_Cont_csd_g3') & \\
				tex \hline
				tex N: & `NPS_Dummy_coef_csdidga_N' & & `DAP_Dummy_coef_csdidga_N' & & `F_Dummy_coef_csdidga_N' & & `NPS_Cont_coef_csdidga_N' & & `DAP_Cont_coef_csdidga_N' & & `F_Cont_coef_csdidga_N' & \\
				tex \hline \hline
				tex \end{tabular}
				tex {\footnotesize \onehalfspacing Notes: The table reports the doubly robust difference in difference estimate of the effect of FBFs on each column's outcome variable. The overall combined effect across cohorts is reported on top, followed by the effect for Cohort 1a (households first treated with a FBF, using not yet treated households as the control and comparing the first and second periods). The Cohort 1b effect is for households first treated with a FBF, using never treated households as the control and comparing the first and third periods. The Cohort 1 combined effect combines the Cohort 1a and 1b effects using the Callaway and Sant'Anna estimator. The Cohort 2 effect reports the effect for households in the 2nd treated FBF group, using never treated households as the control and comparing the second and third periods. For each coefficient of interest, we report our preferred standard errors first in parentheses (wild bootstrapped standard errors clustered at the zone level), followed by wild bootstrapped standard errors clustered at the \textit{woreda} level in brackets. Stars reported alongside the standard errors represent the statistical significance of the corresponding coefficient: $* \left(p<0.1\right), ** \left(p<0.05\right), *** \left(p<0.01\right)$.}
				tex \end{threeparttable}
				tex \end{adjustbox}
				tex \end{table}
				texdoc close 


	//Table S.3: Effect of FBFs on Yields and DAP prices overall and by cohort - Doubly Robust Difference in Difference Model
				cap texdoc close
				texdoc init "Tables/YieldPrice_main_csdid", replace force 
				tex \begin{table}
				tex \begin{adjustbox}{width=.7\linewidth, keepaspectratio}
				tex \begin{threeparttable}
				tex \caption{Effect of FBFs on Yields and DAP prices overall and by cohort - Doubly Robust Difference in Difference Model}\label{results_csdid_YieldPrice}
				tex \begin{tabular}{lclclclcl/*cl*/}
				tex \hline \hline
				tex & \multicolumn{2}{c}{(1)} & \multicolumn{2}{c}{(2)} & \multicolumn{2}{c}{(3)} & \multicolumn{2}{c}{(4)} /*& \multicolumn{2}{c}{(5)}*/ \\
				tex & \multicolumn{2}{c}{Teff} & \multicolumn{2}{c}{Wheat} & \multicolumn{2}{c}{Maize} /*& \multicolumn{2}{c}{Area}*/ & \multicolumn{2}{c}{DAP} \\
				tex & \multicolumn{2}{c}{yields} & \multicolumn{2}{c}{yields} & \multicolumn{2}{c}{yields} /*& \multicolumn{2}{c}{harvested}*/ & \multicolumn{2}{c}{price} \\
				tex & \multicolumn{2}{c}{(kg/ha)} & \multicolumn{2}{c}{(kg/ha)} & \multicolumn{2}{c}{(kg/ha)} /*& \multicolumn{2}{c}{(ha)}*/ & \multicolumn{2}{c}{(ETB/kg)} \\
				tex \hline
				tex Overall & `Yield_1_coef_csdidga' &  & `Yield_2_coef_csdidga' & & `Yield_3_coef_csdidga' /*& & `Area_t_coef_csdidga'*/ & & `DAP_pr_coef_csdidga' & \\
				tex & (`Yield_1_coef_csdidga_se_A05') &  `Yield_1_coef_csdidga_st_A05' & (`Yield_2_coef_csdidga_se_A05') &  `Yield_2_coef_csdidga_st_A05' & (`Yield_3_coef_csdidga_se_A05') &  `Yield_3_coef_csdidga_st_A05' /*& (`Area_t_coef_csdidga_se_A05') &  `Area_t_coef_csdidga_st_A05'*/ & (`DAP_pr_coef_csdidga_se_A05') &  `DAP_pr_coef_csdidga_st_A05' \\
				tex & [`Yield_1_coef_csdidga_se_A04'] &  `Yield_1_coef_csdidga_st_A04' & [`Yield_2_coef_csdidga_se_A04'] &  `Yield_2_coef_csdidga_st_A04' & [`Yield_3_coef_csdidga_se_A04'] &  `Yield_3_coef_csdidga_st_A04' /*& [`Area_t_coef_csdidga_se_A04'] &  `Area_t_coef_csdidga_st_A04'*/ & [`DAP_pr_coef_csdidga_se_A04'] &  `DAP_pr_coef_csdidga_st_A04' \\
				tex Cohort 1a & `Yield_1_coef_csdidc1a' &  & `Yield_2_coef_csdidc1a' & & `Yield_3_coef_csdidc1a' /*& & `Area_t_coef_csdidc1a'*/ & & `DAP_pr_coef_csdidc1a' & \\
				tex & (`Yield_1_coef_csdidc1a_se_A05') &  `Yield_1_coef_csdidc1a_st_A05' & (`Yield_2_coef_csdidc1a_se_A05') &  `Yield_2_coef_csdidc1a_st_A05' & (`Yield_3_coef_csdidc1a_se_A05') &  `Yield_3_coef_csdidc1a_st_A05' /*& (`Area_t_coef_csdidc1a_se_A05') &  `Area_t_coef_csdidc1a_st_A05'*/ & (`DAP_pr_coef_csdidc1a_se_A05') &  `DAP_pr_coef_csdidc1a_st_A05' \\
				tex & [`Yield_1_coef_csdidc1a_se_A04'] &  `Yield_1_coef_csdidc1a_st_A04' & [`Yield_2_coef_csdidc1a_se_A04'] &  `Yield_2_coef_csdidc1a_st_A04' & [`Yield_3_coef_csdidc1a_se_A04'] &  `Yield_3_coef_csdidc1a_st_A04' /*& [`Area_t_coef_csdidc1a_se_A04'] &  `Area_t_coef_csdidc1a_st_A04'*/ & [`DAP_pr_coef_csdidc1a_se_A04'] &  `DAP_pr_coef_csdidc1a_st_A04' \\
				tex Cohort 1b & `Yield_1_coef_csdidc1b' &  & `Yield_2_coef_csdidc1b' & & `Yield_3_coef_csdidc1b' /*& & `Area_t_coef_csdidc1b'*/ & & `DAP_pr_coef_csdidc1b'  & \\
				tex & (`Yield_1_coef_csdidc1b_se_A05') &  `Yield_1_coef_csdidc1b_st_A05' & (`Yield_2_coef_csdidc1b_se_A05') &  `Yield_2_coef_csdidc1b_st_A05' & (`Yield_3_coef_csdidc1b_se_A05') &  `Yield_3_coef_csdidc1b_st_A05' /*& (`Area_t_coef_csdidc1b_se_A05') &  `Area_t_coef_csdidc1b_st_A05'*/ & (`DAP_pr_coef_csdidc1b_se_A05') &  `DAP_pr_coef_csdidc1b_st_A05' \\
				tex & [`Yield_1_coef_csdidc1b_se_A04'] &  `Yield_1_coef_csdidc1b_st_A04' & [`Yield_2_coef_csdidc1b_se_A04'] &  `Yield_2_coef_csdidc1b_st_A04' & [`Yield_3_coef_csdidc1b_se_A04'] &  `Yield_3_coef_csdidc1b_st_A04' /*& [`Area_t_coef_csdidc1b_se_A04'] &  `Area_t_coef_csdidc1b_st_A04'*/ & [`DAP_pr_coef_csdidc1b_se_A04'] &  `DAP_pr_coef_csdidc1b_st_A04' \\
				tex Cohort 1 combined & `Yield_1_coef_csdidg2' &  & `Yield_2_coef_csdidg2' & & `Yield_3_coef_csdidg2' /*& & `Area_t_coef_csdidg2'*/ & & `DAP_pr_coef_csdidg2'  & \\
				tex & (`Yield_1_coef_csdidg2_se_A05') &  `Yield_1_coef_csdidg2_st_A05' & (`Yield_2_coef_csdidg2_se_A05') &  `Yield_2_coef_csdidg2_st_A05' & (`Yield_3_coef_csdidg2_se_A05') &  `Yield_3_coef_csdidg2_st_A05' /*& (`Area_t_coef_csdidg2_se_A05') &  `Area_t_coef_csdidg2_st_A05'*/ & (`DAP_pr_coef_csdidg2_se_A05') &  `DAP_pr_coef_csdidg2_st_A05' \\
				tex & [`Yield_1_coef_csdidg2_se_A04'] &  `Yield_1_coef_csdidg2_st_A04' & [`Yield_2_coef_csdidg2_se_A04'] &  `Yield_2_coef_csdidg2_st_A04' & [`Yield_3_coef_csdidg2_se_A04'] &  `Yield_3_coef_csdidg2_st_A04' /*& [`Area_t_coef_csdidg2_se_A04'] &  `Area_t_coef_csdidg2_st_A04'*/ & [`DAP_pr_coef_csdidg2_se_A04'] &  `DAP_pr_coef_csdidg2_st_A04' \\
				tex Cohort 2 & `Yield_1_coef_csdidg3' &  & `Yield_2_coef_csdidg3' & & `Yield_3_coef_csdidg3' /*& & `Area_t_coef_csdidg3'*/ & & `DAP_pr_coef_csdidg3' & \\
				tex & (`Yield_1_coef_csdidg3_se_A05') &  `Yield_1_coef_csdidg3_st_A05' & (`Yield_2_coef_csdidg3_se_A05') &  `Yield_2_coef_csdidg3_st_A05' & (`Yield_3_coef_csdidg3_se_A05') &  `Yield_3_coef_csdidg3_st_A05' /*& (`Area_t_coef_csdidg3_se_A05') &  `Area_t_coef_csdidg3_st_A05'*/ & (`DAP_pr_coef_csdidg3_se_A05') &  `DAP_pr_coef_csdidg3_st_A05' \\
				tex & [`Yield_1_coef_csdidg3_se_A04'] &  `Yield_1_coef_csdidg3_st_A04' & [`Yield_2_coef_csdidg3_se_A04'] &  `Yield_2_coef_csdidg3_st_A04' & [`Yield_3_coef_csdidg3_se_A04'] &  `Yield_3_coef_csdidg3_st_A04' /*& [`Area_t_coef_csdidg3_se_A04'] &  `Area_t_coef_csdidg3_st_A04'*/ & [`DAP_pr_coef_csdidg3_se_A04'] &  `DAP_pr_coef_csdidg3_st_A04' \\ 
				tex \hline 
				tex Pre-period control mean (sd) & & & & & & & & /*& &*/  \\
				tex \hspace{3mm} Overall & `Yield_1_cmean_ga' & & `Yield_2_cmean_ga' & & `Yield_3_cmean_ga' /*& & `Area_t_cmean_ga'*/ & & `DAP_pr_cmean_ga' & \\
				tex  & (`Yield_1_csd_ga') & & (`Yield_2_csd_ga') & & (`Yield_3_csd_ga') /*& & (`Area_t_csd_ga')*/ & & (`DAP_pr_csd_ga') & \\
				tex \hspace{3mm} Cohort 1a & `Yield_1_cmean_c1a' & & `Yield_2_cmean_c1a' & & `Yield_3_cmean_c1a' /*& & `Area_t_cmean_c1a'*/ & & `DAP_pr_cmean_c1a' & \\
				tex & (`Yield_1_csd_c1a') & & (`Yield_2_csd_c1a') & & (`Yield_3_csd_c1a') /*& & (`Area_t_csd_c1a')*/ & & (`DAP_pr_csd_c1a')  & \\
				tex \hspace{3mm} Cohort 1b & `Yield_1_cmean_c1b' & & `Yield_2_cmean_c1b' & & `Yield_3_cmean_c1b' /*& & `Area_t_cmean_c1b'*/ & & `DAP_pr_cmean_c1b'  & \\
				tex & (`Yield_1_csd_c1b') & & (`Yield_2_csd_c1b') & & (`Yield_3_csd_c1b') /*& & (`Area_t_csd_c1b')*/ & & (`DAP_pr_csd_c1b') & \\
				tex \hspace{3mm} Cohort 1 combined & `Yield_1_cmean_g2' & & `Yield_2_cmean_g2' & & `Yield_3_cmean_g2' /*& & `Area_t_cmean_g2'*/ & & `DAP_pr_cmean_g2' & \\
				tex  & (`Yield_1_csd_g2') & & (`Yield_2_csd_g2') & & (`Yield_3_csd_g2') /*& & (`Area_t_csd_g2')*/ & & (`DAP_pr_csd_g2') & \\
				tex \hspace{3mm} Cohort 2 & `Yield_1_cmean_g3' & & `Yield_2_cmean_g3' & & `Yield_3_cmean_g3' /*& & `Area_t_cmean_g3'*/ & & `DAP_pr_cmean_g3' & \\
				tex  & (`Yield_1_csd_g3') & & (`Yield_2_csd_g3') & & (`Yield_3_csd_g3') /*& & (`Area_t_csd_g3')*/ & & (`DAP_pr_csd_g3') & \\
				tex \hline
				tex N: & `Yield_1_coef_csdidga_N' & & `Yield_2_coef_csdidga_N' & & `Yield_3_coef_csdidga_N' /*& & `Area_t_coef_csdidga_N'*/ & & `DAP_pr_coef_csdidga_N' & \\
				tex \hline \hline
				tex \end{tabular}
				tex {\footnotesize \onehalfspacing Notes: The table reports the doubly robust difference in difference estimate of the effect of FBFs on each column's outcome variable. The overall combined effect across cohorts is reported on top, followed by the effect for Cohort 1a (households first treated with a FBF, using not yet treated households as the control and comparing the first and second periods). The Cohort 1b effect is for households first treated with a FBF, using never treated households as the control and comparing the first and third periods. The Cohort 1 combined effect combines the Cohort 1a and 1b effects using the Callaway and Sant'Anna estimator. The Cohort 2 effect reports the effect for households in the 2nd treated FBF group, using never treated households as the control and comparing the second and third periods. For each coefficient of interest, we report our preferred standard errors first in parentheses (wild bootstrapped standard errors clustered at the zone level), followed by wild bootstrapped standard errors clustered at the \textit{woreda} level in brackets. Stars reported alongside the standard errors represent the statistical significance of the corresponding coefficient: $* \left(p<0.1\right), ** \left(p<0.05\right), *** \left(p<0.01\right)$.}
				tex \end{threeparttable}
				tex \end{adjustbox}
				tex \end{table}
				texdoc close 


			//Table S.4: Effect of FBFs on Gross Crop Revenue and Household Consumption overall and by cohort - Doubly Robust Difference in Difference Model

				cap texdoc close
				texdoc init "Tables/ValueHarvestCons_main_csdid", replace force 
				tex \begin{table}
				tex \begin{adjustbox}{width=.7\linewidth, keepaspectratio}
				tex \begin{threeparttable}
				tex \caption{Effect of FBFs on Gross Crop Revenue and Household Consumption overall and by cohort - Doubly Robust Difference in Difference Model}\label{results_csdid_ValueHarvestCons}
				tex \begin{tabular}{lclclclclcl}
				tex \hline \hline
				tex & \multicolumn{2}{c}{(1)} & \multicolumn{2}{c}{(2)} & \multicolumn{2}{c}{(3)} & \multicolumn{2}{c}{(4)} & \multicolumn{2}{c}{(5)} \\
				tex & \multicolumn{2}{c}{Total} & \multicolumn{2}{c}{Teff} & \multicolumn{2}{c}{Wheat} & \multicolumn{2}{c}{Maize} & \multicolumn{2}{c}{Consumption} \\
				tex & \multicolumn{2}{c}{value} & \multicolumn{2}{c}{value} & \multicolumn{2}{c}{value} & \multicolumn{2}{c}{value} & \multicolumn{2}{c}{per adult} \\
				tex & \multicolumn{2}{c}{harvested} & \multicolumn{2}{c}{harvested} & \multicolumn{2}{c}{harvested} & \multicolumn{2}{c}{harvested} & \multicolumn{2}{c}{equivalent} \\
				tex & \multicolumn{2}{c}{(ETB)} & \multicolumn{2}{c}{(ETB)} & \multicolumn{2}{c}{(ETB)} & \multicolumn{2}{c}{(ETB)} & \multicolumn{2}{c}{(2011 ETB/day)} \\
				tex \hline
				tex Overall & `v_harv_coef_csdidga' &  & `v_harv1_coef_csdidga' & & `v_harv2_coef_csdidga' & & `v_harv3_coef_csdidga' & & `pcexpaeDR_coef_csdidga' & \\
				tex & (`v_harv_coef_csdidga_se_A05') &  `v_harv_coef_csdidga_st_A05' & (`v_harv1_coef_csdidga_se_A05') &  `v_harv1_coef_csdidga_st_A05' & (`v_harv2_coef_csdidga_se_A05') &  `v_harv2_coef_csdidga_st_A05' & (`v_harv3_coef_csdidga_se_A05') &  `v_harv3_coef_csdidga_st_A05' & (`pcexpaeDR_coef_csdidga_se_A05') &  `pcexpaeDR_coef_csdidga_st_A05' \\
				tex & [`v_harv_coef_csdidga_se_A04'] &  `v_harv_coef_csdidga_st_A04' & [`v_harv1_coef_csdidga_se_A04'] &  `v_harv1_coef_csdidga_st_A04' & [`v_harv2_coef_csdidga_se_A04'] &  `v_harv2_coef_csdidga_st_A04' & [`v_harv3_coef_csdidga_se_A04'] &  `v_harv3_coef_csdidga_st_A04' & [`pcexpaeDR_coef_csdidga_se_A04'] &  `pcexpaeDR_coef_csdidga_st_A04' \\
				tex Cohort 1a & `v_harv_coef_csdidc1a' &  & `v_harv1_coef_csdidc1a' & & `v_harv2_coef_csdidc1a' & & `v_harv3_coef_csdidc1a' & & `pcexpaeDR_coef_csdidc1a' & \\
				tex & (`v_harv_coef_csdidc1a_se_A05') &  `v_harv_coef_csdidc1a_st_A05' & (`v_harv1_coef_csdidc1a_se_A05') &  `v_harv1_coef_csdidc1a_st_A05' & (`v_harv2_coef_csdidc1a_se_A05') &  `v_harv2_coef_csdidc1a_st_A05' & (`v_harv3_coef_csdidc1a_se_A05') &  `v_harv3_coef_csdidc1a_st_A05' & (`pcexpaeDR_coef_csdidc1a_se_A05') &  `pcexpaeDR_coef_csdidc1a_st_A05' \\
				tex & [`v_harv_coef_csdidc1a_se_A04'] &  `v_harv_coef_csdidc1a_st_A04' & [`v_harv1_coef_csdidc1a_se_A04'] &  `v_harv1_coef_csdidc1a_st_A04' & [`v_harv2_coef_csdidc1a_se_A04'] &  `v_harv2_coef_csdidc1a_st_A04' & [`v_harv3_coef_csdidc1a_se_A04'] &  `v_harv3_coef_csdidc1a_st_A04' & [`pcexpaeDR_coef_csdidc1a_se_A04'] &  `pcexpaeDR_coef_csdidc1a_st_A04' \\
				tex Cohort 1b & `v_harv_coef_csdidc1b' &  & `v_harv1_coef_csdidc1b' & & `v_harv2_coef_csdidc1b' & & `v_harv3_coef_csdidc1b' & & `pcexpaeDR_coef_csdidc1b'  & \\
				tex & (`v_harv_coef_csdidc1b_se_A05') &  `v_harv_coef_csdidc1b_st_A05' & (`v_harv1_coef_csdidc1b_se_A05') &  `v_harv1_coef_csdidc1b_st_A05' & (`v_harv2_coef_csdidc1b_se_A05') &  `v_harv2_coef_csdidc1b_st_A05' & (`v_harv3_coef_csdidc1b_se_A05') &  `v_harv3_coef_csdidc1b_st_A05' & (`pcexpaeDR_coef_csdidc1b_se_A05') &  `pcexpaeDR_coef_csdidc1b_st_A05' \\
				tex & [`v_harv_coef_csdidc1b_se_A04'] &  `v_harv_coef_csdidc1b_st_A04' & [`v_harv1_coef_csdidc1b_se_A04'] &  `v_harv1_coef_csdidc1b_st_A04' & [`v_harv2_coef_csdidc1b_se_A04'] &  `v_harv2_coef_csdidc1b_st_A04' & [`v_harv3_coef_csdidc1b_se_A04'] &  `v_harv3_coef_csdidc1b_st_A04' & [`pcexpaeDR_coef_csdidc1b_se_A04'] &  `pcexpaeDR_coef_csdidc1b_st_A04' \\
				tex Cohort 1 combined & `v_harv_coef_csdidg2' &  & `v_harv1_coef_csdidg2' & & `v_harv2_coef_csdidg2' & & `v_harv3_coef_csdidg2' & & `pcexpaeDR_coef_csdidg2'  & \\
				tex & (`v_harv_coef_csdidg2_se_A05') &  `v_harv_coef_csdidg2_st_A05' & (`v_harv1_coef_csdidg2_se_A05') &  `v_harv1_coef_csdidg2_st_A05' & (`v_harv2_coef_csdidg2_se_A05') &  `v_harv2_coef_csdidg2_st_A05' & (`v_harv3_coef_csdidg2_se_A05') &  `v_harv3_coef_csdidg2_st_A05' & (`pcexpaeDR_coef_csdidg2_se_A05') &  `pcexpaeDR_coef_csdidg2_st_A05' \\
				tex & [`v_harv_coef_csdidg2_se_A04'] &  `v_harv_coef_csdidg2_st_A04' & [`v_harv1_coef_csdidg2_se_A04'] &  `v_harv1_coef_csdidg2_st_A04' & [`v_harv2_coef_csdidg2_se_A04'] &  `v_harv2_coef_csdidg2_st_A04' & [`v_harv3_coef_csdidg2_se_A04'] &  `v_harv3_coef_csdidg2_st_A04' & [`pcexpaeDR_coef_csdidg2_se_A04'] &  `pcexpaeDR_coef_csdidg2_st_A04' \\
				tex Cohort 2 & `v_harv_coef_csdidg3' &  & `v_harv1_coef_csdidg3' & & `v_harv2_coef_csdidg3' & & `v_harv3_coef_csdidg3' & & `pcexpaeDR_coef_csdidg3' & \\
				tex & (`v_harv_coef_csdidg3_se_A05') &  `v_harv_coef_csdidg3_st_A05' & (`v_harv1_coef_csdidg3_se_A05') &  `v_harv1_coef_csdidg3_st_A05' & (`v_harv2_coef_csdidg3_se_A05') &  `v_harv2_coef_csdidg3_st_A05' & (`v_harv3_coef_csdidg3_se_A05') &  `v_harv3_coef_csdidg3_st_A05' & (`pcexpaeDR_coef_csdidg3_se_A05') &  `pcexpaeDR_coef_csdidg3_st_A05' \\
				tex & [`v_harv_coef_csdidg3_se_A04'] &  `v_harv_coef_csdidg3_st_A04' & [`v_harv1_coef_csdidg3_se_A04'] &  `v_harv1_coef_csdidg3_st_A04' & [`v_harv2_coef_csdidg3_se_A04'] &  `v_harv2_coef_csdidg3_st_A04' & [`v_harv3_coef_csdidg3_se_A04'] &  `v_harv3_coef_csdidg3_st_A04' & [`pcexpaeDR_coef_csdidg3_se_A04'] &  `pcexpaeDR_coef_csdidg3_st_A04' \\ 
				tex \hline 
				tex Pre-period control mean (sd) & & & & & & & & & &  \\
				tex \hspace{3mm} Overall & `v_harv_cmean_ga' & & `v_harv1_cmean_ga' & & `v_harv2_cmean_ga' & & `v_harv3_cmean_ga' & & `pcexpaeDR_cmean_ga' & \\
				tex  & (`v_harv_csd_ga') & & (`v_harv1_csd_ga') & & (`v_harv2_csd_ga') & & (`v_harv3_csd_ga') & & (`pcexpaeDR_csd_ga') & \\
				tex \hspace{3mm} Cohort 1a & `v_harv_cmean_c1a' & & `v_harv1_cmean_c1a' & & `v_harv2_cmean_c1a' & & `v_harv3_cmean_c1a' & & `pcexpaeDR_cmean_c1a' & \\
				tex & (`v_harv_csd_c1a') & & (`v_harv1_csd_c1a') & & (`v_harv2_csd_c1a') & & (`v_harv3_csd_c1a') & & (`pcexpaeDR_csd_c1a')  & \\
				tex \hspace{3mm} Cohort 1b & `v_harv_cmean_c1b' & & `v_harv1_cmean_c1b' & & `v_harv2_cmean_c1b' & & `v_harv3_cmean_c1b' & & `pcexpaeDR_cmean_c1b'  & \\
				tex & (`v_harv_csd_c1b') & & (`v_harv1_csd_c1b') & & (`v_harv2_csd_c1b') & & (`v_harv3_csd_c1b') & & (`pcexpaeDR_csd_c1b') & \\
				tex \hspace{3mm} Cohort 1 combined & `v_harv_cmean_g2' & & `v_harv1_cmean_g2' & & `v_harv2_cmean_g2' & & `v_harv3_cmean_g2' & & `pcexpaeDR_cmean_g2' & \\
				tex  & (`v_harv_csd_g2') & & (`v_harv1_csd_g2') & & (`v_harv2_csd_g2') & & (`v_harv3_csd_g2') & & (`pcexpaeDR_csd_g2') & \\
				tex \hspace{3mm} Cohort 2 & `v_harv_cmean_g3' & & `v_harv1_cmean_g3' & & `v_harv2_cmean_g3' & & `v_harv3_cmean_g3' & & `pcexpaeDR_cmean_g3' & \\
				tex  & (`v_harv_csd_g3') & & (`v_harv1_csd_g3') & & (`v_harv2_csd_g3') & & (`v_harv3_csd_g3') & & (`pcexpaeDR_csd_g3') & \\
				tex \hline
				tex N: & `v_harv_coef_csdidga_N' & & `v_harv1_coef_csdidga_N' & & `v_harv2_coef_csdidga_N' & & `v_harv3_coef_csdidga_N' & & `pcexpaeDR_coef_csdidga_N' & \\
				tex \hline \hline
				tex \end{tabular}
				tex {\footnotesize \onehalfspacing Notes: The table reports the doubly robust difference in difference estimate of the effect of FBFs on each column's outcome variable. The overall combined effect across cohorts is reported on top, followed by the effect for Cohort 1a (households first treated with a FBF, using not yet treated households as the control and comparing the first and second periods). The Cohort 1b effect is for households first treated with a FBF, using never treated households as the control and comparing the first and third periods. The Cohort 1 combined effect combines the Cohort 1a and 1b effects using the Callaway and Sant'Anna estimator. The Cohort 2 effect reports the effect for households in the 2nd treated FBF group, using never treated households as the control and comparing the second and third periods. For each coefficient of interest, we report our preferred standard errors first in parentheses (wild bootstrapped standard errors clustered at the zone level), followed by wild bootstrapped standard errors clustered at the \textit{woreda} level in brackets. Stars reported alongside the standard errors represent the statistical significance of the corresponding coefficient: $* \left(p<0.1\right), ** \left(p<0.05\right), *** \left(p<0.01\right)$.}
				tex \end{threeparttable}
				tex \end{adjustbox}
				tex \end{table}


			//Table S.5, - S.19: DID Tables for each outcome - outcomelist includes NPS_Dummy (for S.5), DAP_Dummy (for S.6), F_Dummy (for S.7), NPS_Cont (for S.8), DAP_Cont (for S.9), F_Cont (S.10), DAP_pr (for S.11), Yield_1 (for S.12), Yield_2 (for S.13), Yield_3(for S.14), v_harv(for S.15), v_harv1 (for S.16), v_harv2 (for S.17), v_harv3 (for S.18),  pcexpaeDR (for S.19)

				
				foreach outcome in `outcomelist' {
				foreach c in 1a 1b 2 {
					use `dataset_c`c'', clear
					
					su `outcome' if (cohort`c'toinclude==1) & (cohort`c'pre==1) & (cohort`c'control==1)
					loc `outcome'_cmean_c`c' = r(mean)
					loc `outcome'_cmean_c`c' : display %09.3fc ``outcome'_cmean_c`c''
					loc `outcome'_csd_c`c' = r(sd)
					loc `outcome'_csd_c`c' : display %09.3fc ``outcome'_csd_c`c''
					
					*Remove household partial effects
					cap drop `outcome'_id `outcome'_dm
					areg `outcome' if cohort`c'toinclude==1, absorb(id)
					predict `outcome'_id
					g `outcome'_dm = `outcome' - `outcome'_id

					
					** Three regressions
					loc r = 1
					wildbootstrap reg `outcome'_dm cohort`c'treatpost  cohort`c'treat cohort`c'post `covars1' if cohort`c'toinclude==1, cluster(`clustsize_default') rseed(`seed')
					mat `outcome'_b_c`c'_r`r' = e(b)'
					mata st_matrix("`outcome'_se_c`c'_r`r'",sqrt(diagonal(st_matrix("e(V)"))))
					loc `outcome'_N_c`c'_r`r' = e(N)
					loc `outcome'_N_c`c'_r`r' : display %09.0fc ``outcome'_N_c`c'_r`r''
					loc `outcome'_Nc_c`c'_r`r' = e(N_clust)
					loc `outcome'_Nc_c`c'_r`r' : display %09.0fc ``outcome'_Nc_c`c'_r`r''
					estimates store `outcome'_c`c'_r`r'
					
					loc r = 2
					wildbootstrap reg `outcome'_dm cohort`c'treatpost  cohort`c'treat cohort`c'post `covars1' if cohort`c'toinclude==1, cluster(`clustsize_alt') rseed(`seed')
					mat `outcome'_b_c`c'_r`r' = e(b)'
					mata st_matrix("`outcome'_se_c`c'_r`r'",sqrt(diagonal(st_matrix("e(V)"))))
					loc `outcome'_N_c`c'_r`r' = e(N)
					loc `outcome'_N_c`c'_r`r' : display %09.0fc ``outcome'_N_c`c'_r`r''
					loc `outcome'_Nc_c`c'_r`r' = e(N_clust)
					loc `outcome'_Nc_c`c'_r`r' : display %09.0fc ``outcome'_Nc_c`c'_r`r''
					estimates store `outcome'_c`c'_r`r'

					loc r = 3
					sort `clustsize_spatial' id year
					reg `outcome'_dm cohort`c'treatpost  cohort`c'treat cohort`c'post `covars1' if cohort`c'toinclude==1, cluster(`clustsize_spatial')
					loc `outcome'_N_c`c'_r`r' = e(N)
					loc `outcome'_N_c`c'_r`r' : display %09.0fc ``outcome'_N_c`c'_r`r''
					loc `outcome'_Nc_c`c'_r`r' = e(N_clust)
					loc `outcome'_Nc_c`c'_r`r' : display %09.0fc ``outcome'_Nc_c`c'_r`r''
					scpc 
					mat results = e(scpcstats)
					mat `outcome'_b_c`c'_r`r' = results[1...,1]
					mat `outcome'_se_c`c'_r`r' = results[1...,2]
					mat list `outcome'_b_c`c'_r`r'
					mat list  `outcome'_se_c`c'_r`r'
					
					loc J = rowsof(`outcome'_b_c`c'_r`r')
					mat `outcome'_st_c`c'_r0  = J(`J',1,0) /*blank stars matrix for the beta rows*/
					forvalues r=1/3 {
						loc J = rowsof(`outcome'_b_c`c'_r`r')
						mat `outcome'_st_c`c'_r`r' = J(`J',1,0)
						mat list `outcome'_st_c`c'_r`r'
						display "ESTIMATES: outcome `outcome'; cohort `c'; estimates `r'."
						forvalues j=1/`J' {
							display "j: `j'"
							loc star = 0 
							loc beta = `outcome'_b_c`c'_r`r'[`j',1]
							display "beta: `beta'"
							loc se = `outcome'_se_c`c'_r`r'[`j',1]
							display "se: `se'"
							loc N = ``outcome'_Nc_c`c'_r`r''
							display "N: `N'"
							foreach sig in `siglist' {
								cap loc tstat = (abs(`beta'/`se')>invttail(`N',`sig'/2))
								display "sig: `sig'; tstat: `tstat'" 
								if `tstat'==float(1) {
									loc star = `star'+1
									}
								}
							display "stars: `star'"
							mat `outcome'_st_c`c'_r`r'[`j',1] = `star'
							}
						}
						
					mat `outcome'_statmat_c`c' = `outcome'_b_c`c'_r1,`outcome'_se_c`c'_r1,`outcome'_se_c`c'_r2,`outcome'_se_c`c'_r3
					mat `outcome'_statmat_c`c'_st = `outcome'_st_c`c'_r0, `outcome'_st_c`c'_r1, `outcome'_st_c`c'_r2, `outcome'_st_c`c'_r3
					}
					
				mat `outcome'_statmat = `outcome'_statmat_c1a, `outcome'_statmat_c1b, `outcome'_statmat_c2
				mat `outcome'_statmat_st = `outcome'_statmat_c1a_st, `outcome'_statmat_c1b_st, `outcome'_statmat_c2_st

				frmttable using "Tables/`outcome'_cohortdid", replace statmat(`outcome'_statmat) sdec(3) /*title("Effect of FBFs on ``outcome'_label': cohort-specific difference in difference outcome regressions")*/ substat(3) ctitle("","Cohort 1a","Cohort 1b", "Cohort 2") varlabels  annotate(`outcome'_statmat_st) asymbol(*,**,***)  tex fragment addrows("N:", "``outcome'_N_c1a_r1'", "``outcome'_N_c1b_r1'", "``outcome'_N_c2_r1'" \ "Pre-period control mean", "``outcome'_cmean_c1a'", "``outcome'_cmean_c1b'", "``outcome'_cmean_c2'" \ "\hspace{3mm} (sd)", "``outcome'_csd_c1a'", "``outcome'_csd_c1b'", "``outcome'_csd_c2'") hlines(11{0}1001)
				}



			//Table S.21: Effect of FBFs on Fertilizer Adoption Overall and by cohort - spillovers analysis
				* identify the near-treat for year 2 and year 3 based on contuous travel time (donut)
				use `dataset_all', clear

				foreach t in 1 2 {
				g treat_`t'_yes = Treatment_`t'_continous < 4 
				g treat_`t'_near = Treatment_`t'_continous >= 4 & Treatment_`t'_continous < 5.5
				g treat_`t'_no = Treatment_`t'_continous >= 5.5
				}

				tab treat_1_yes year
				tab treat_1_near year
				tab treat_1_no year

				tempfile dataset_interim_neartreat
				save `dataset_interim_neartreat', replace


				* identify the near-treat for year 2 and year 3 based on woredas with a mixture 

				use `dataset_all', clear

				bysort `clustsize_alt': egen Treatment_1_`clustsize_alt' = mean(Treatment_1)
				bysort `clustsize_alt': egen Treatment_2_`clustsize_alt' = mean(Treatment_2)
				tab Treatment_1 if cohort1atoinclude==1
				tab Treatment_1 if cohort1btoinclude==1
				tab Treatment_2 if cohort2toinclude==1

				foreach t in 1 2 {
				g treat_`t'_yes = Treatment_`t'_`clustsize_alt'== 1  
				g treat_`t'_near = Treatment_`t'_`clustsize_alt' > 0  & Treatment_`t'_`clustsize_alt' < 1 /*1680 sample for T1 and 1773 for T2, a few are barely above zero, a few barely below 1, some in the middle */
				g treat_`t'_no = Treatment_`t'_`clustsize_alt' == 0
				}

				tab treat_1_yes year
				tab treat_1_near year
				tab treat_1_no year

				tempfile dataset_interim_A04
				save `dataset_interim_A04', replace


				** First Cohort periods 1 & 2 DID setup -- excluding near treat
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_1_near==1
				g nnt_cohort1atreat = treat_1_yes==1
				g nnt_cohort1acontrol = treat_1_no==1 
				g nnt_cohort1apre = (year==1)
				g nnt_cohort1apost = (year==2)
				g nnt_cohort1atreatpost = nnt_cohort1atreat * nnt_cohort1apost
				g nnt_cohort1atoinclude = (nnt_cohort1atreat==1 | nnt_cohort1acontrol==1) & (nnt_cohort1apre==1 | nnt_cohort1apost==1)
				keep if nnt_cohort1atoinclude==1
				tempfile nnt_cohort1a
				save `nnt_cohort1a', replace


				** First Cohort periods 1 & 2 DID setup -- spillovers
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_1_yes==1
				g nto_cohort1atreat = treat_1_near==1
				g nto_cohort1acontrol = treat_1_no==1 
				g nto_cohort1apre = (year==1)
				g nto_cohort1apost = (year==2)
				g nto_cohort1atreatpost = nto_cohort1atreat * nto_cohort1apost
				g nto_cohort1atoinclude = (nto_cohort1atreat==1 | nto_cohort1acontrol==1) & (nto_cohort1apre==1 | nto_cohort1apost==1)
				keep if nto_cohort1atoinclude==1
				tempfile nto_cohort1a
				save `nto_cohort1a', replace


				** First Cohort periods 1 & 3 DID setup - excluding near treat
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_1_near==1
				g nnt_cohort1btreat = treat_1_yes==1
				g nnt_cohort1bcontrol = treat_1_no==1 & treat_2_no==1 /*Possibly could include treat_2_near?*/
				g nnt_cohort1bpre = (year==1)
				g nnt_cohort1bpost = (year==3)
				g nnt_cohort1btreatpost = nnt_cohort1btreat * nnt_cohort1bpost
				g nnt_cohort1btoinclude = (nnt_cohort1btreat==1 | nnt_cohort1bcontrol==1) & (nnt_cohort1bpre==1 | nnt_cohort1bpost==1)
				keep if nnt_cohort1btoinclude==1
				tempfile nnt_cohort1b
				save `nnt_cohort1b', replace


				** First Cohort periods 1 & 3 DID setup - spillovers
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_1_yes==1
				g nto_cohort1btreat = treat_1_near==1
				g nto_cohort1bcontrol = treat_1_no==1 & treat_2_no==1 /*Possibly could include treat_2_near?*/
				g nto_cohort1bpre = (year==1)
				g nto_cohort1bpost = (year==3)
				g nto_cohort1btreatpost = nto_cohort1btreat * nto_cohort1bpost
				g nto_cohort1btoinclude = (nto_cohort1btreat==1 | nto_cohort1bcontrol==1) & (nto_cohort1bpre==1 | nto_cohort1bpost==1)
				keep if nto_cohort1btoinclude==1
				tempfile nto_cohort1b
				save `nto_cohort1b', replace


				** Second period DID setup - excluding near treat
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_2_near==1
				g nnt_cohort2treat = (treat_2_yes==1)
				g nnt_cohort2control = (treat_2_no==1) & treat_1_no==1 /*Possibly could include treat_1_near?*/
				g nnt_cohort2pre = (year==2)
				g nnt_cohort2post = (year==3)
				g nnt_cohort2treatpost = nnt_cohort2treat * nnt_cohort2post
				g nnt_cohort2toinclude = (nnt_cohort2treat==1 | nnt_cohort2control==1) & (nnt_cohort2pre==1 | nnt_cohort2post==1)
				keep if nnt_cohort2toinclude==1
				tempfile nnt_cohort2
				save `nnt_cohort2', replace


				** Second period DID setup - spillovers
				*use `dataset_interim_neartreat', clear
				use `dataset_interim_A04', clear
				drop if treat_2_yes==1
				g nto_cohort2treat = (treat_2_near==1)
				g nto_cohort2control = (treat_2_no==1) & treat_1_no==1 /*Possibly could include treat_1_near?*/
				g nto_cohort2pre = (year==2)
				g nto_cohort2post = (year==3)
				g nto_cohort2treatpost = nto_cohort2treat * nto_cohort2post
				g nto_cohort2toinclude = (nto_cohort2treat==1 | nto_cohort2control==1) & (nto_cohort2pre==1 | nto_cohort2post==1)
				keep if nto_cohort2toinclude==1
				tempfile nto_cohort2
				save `nto_cohort2', replace


				cap texdoc close
				texdoc init "Tables/FertDummyCont_spillovers_drdid", replace force 
				tex \begin{table}
				tex \begin{adjustbox}{width=\linewidth, keepaspectratio}
				tex \begin{threeparttable}
				tex \caption{Effect of FBFs on Fertilizer Adoption Overall and by cohort - spillovers analysis}\label{spillovers_drdid_FertDummyCont}
				tex \begin{tabular}{llcccccccccccc}
				tex \hline \hline
				tex & & \multicolumn{6}{c}{Excluding Near Treated Households} & \multicolumn{6}{c}{Near Treated Households Only} \\
				tex & & \multicolumn{2}{c}{Cohort} & \multicolumn{2}{c}{Cohort} & \multicolumn{2}{c}{Cohort} & \multicolumn{2}{c}{Cohort} & \multicolumn{2}{c}{Cohort} & \multicolumn{2}{c}{Cohort} \\
				tex & & \multicolumn{2}{c}{1a} & \multicolumn{2}{c}{1b} & \multicolumn{2}{c}{2} & \multicolumn{2}{c}{1a} & \multicolumn{2}{c}{1b} & \multicolumn{2}{c}{2} \\
				tex & & \multicolumn{2}{c}{(1)} & \multicolumn{2}{c}{(2)} & \multicolumn{2}{c}{(3)} & \multicolumn{2}{c}{(4)} & \multicolumn{2}{c}{(5)} & \multicolumn{2}{c}{(6)} \\

				tex \hline

				** DRDID results for main tables but with spillovers
				foreach outcome in `outcomelist' {
				foreach c in 1a 1b 2 {
					
					** Two samples -- first excluding near treated hhs then next for only near treated
					foreach sam in nnt nto {
						use ``sam'_cohort`c'', clear
						
						su `outcome' if (`sam'_cohort`c'toinclude==1) & (`sam'_cohort`c'pre==1) & (`sam'_cohort`c'control==1)
						loc `outcome'_`sam'_cmean_c`c' = r(mean)
						loc `outcome'_`sam'_cmean_c`c' : display %09.3fc ``outcome'_`sam'_cmean_c`c''
						loc `outcome'_`sam'_csd_c`c' = r(sd)
						loc `outcome'_`sam'_csd_c`c' : display %09.3fc ``outcome'_`sam'_csd_c`c''
						
						*Remove household partial effects
						cap drop `outcome'_id `outcome'_dm
						areg `outcome', absorb(id)
						predict `outcome'_id
						g `outcome'_dm = `outcome' - `outcome'_id

						
						** Cohort regression
						drdid `outcome'_dm `covars1' if `sam'_cohort`c'toinclude==1, dripw  wboot cluster(`clustsize_default') rseed(`seed') /*ivar(id)*/ time(year) treatment(`sam'_cohort`c'treat) level(95) stub(`outcome'_nnt_c`c'_) noisily 
						*wildbootstrap reg `outcome'_dm cohort`c'treatpost  cohort`c'treat cohort`c'post `covars1' if cohort`c'toinclude==1, cluster(`clustsize_default') rseed(`seed')
						mat beta = e(b)
						mat var = e(V)
						loc `outcome'_`sam'_b_c`c' = beta[1,1]
						loc `outcome'_`sam'_b_c`c' : display %09.3fc ``outcome'_`sam'_b_c`c''
						loc `outcome'_`sam'_se_c`c' = sqrt(var[1,1])
						loc `outcome'_`sam'_se_c`c' : display %09.3fc ``outcome'_`sam'_se_c`c''
						loc `outcome'_`sam'_N_c`c' = e(N)
						loc `outcome'_`sam'_N_c`c' : display %09.0fc ``outcome'_`sam'_N_c`c''
						loc `outcome'_`sam'_Nc_c`c' = e(N_clust)
						loc `outcome'_`sam'_Nc_c`c' : display %09.0fc ``outcome'_`sam'_Nc_c`c''
						estimates store `outcome'_nnt_c`c'
						
						display "ESTIMATES: outcome `outcome'; cohort `c';"

						loc star ""
						cap loc beta = ``outcome'_`sam'_b_c`c''
						display "beta: `beta'"
						cap loc se = ``outcome'_`sam'_se_c`c''
						display "se: `se'"
						cap loc N = ``outcome'_`sam'_Nc_c`c''
						display "N: `N'"
						foreach sig in `siglist' {
							cap loc tstat = (abs(`beta'/`se')>invttail(`N',`sig'/2))
							display "sig: `sig'; tstat: `tstat'" 
							if `tstat'==float(1) {
								loc star `star'*
								}
							}
						display "stars: `star'"
						loc `outcome'_`sam'_st_c`c' `star'
						}
							
					}
				tex ``outcome'_Label' & Treat X Post & ``outcome'_nnt_b_c1a' & ``outcome'_nnt_st_c1a' & ``outcome'_nnt_b_c1b' & ``outcome'_nnt_st_c1b' & ``outcome'_nnt_b_c2' & ``outcome'_nnt_st_c2' & ``outcome'_nto_b_c1a' & ``outcome'_nto_st_c1a' & ``outcome'_nto_b_c1b' & ``outcome'_nto_st_c1b' & ``outcome'_nto_b_c2' & ``outcome'_nto_st_c2' \\
				tex  &  & (``outcome'_nnt_se_c1a') & & (``outcome'_nnt_se_c1b') &  & (``outcome'_nnt_se_c2') &  & (``outcome'_nto_se_c1a') &  & (``outcome'_nto_se_c1b') &  & (``outcome'_nto_se_c2') &  \\
				tex  & Control mean & ``outcome'_nnt_cmean_c1a' & & ``outcome'_nnt_cmean_c1b' &  & ``outcome'_nnt_cmean_c2' &  & ``outcome'_nto_cmean_c1a' &  & ``outcome'_nto_cmean_c1b' &  & ``outcome'_nto_cmean_c2' & \\
				tex  &  & (``outcome'_nnt_csd_c1a') & & (``outcome'_nnt_csd_c1b') &  & (``outcome'_nnt_csd_c2') &  & (``outcome'_nto_csd_c1a') &  & (``outcome'_nto_csd_c1b') &  & (``outcome'_nto_csd_c2') &  \\

				}

				tex \hline
				tex N & & `NPS_Dummy_nnt_N_c1a' & &  `NPS_Dummy_nnt_N_c1b' & & `NPS_Dummy_nnt_N_c2' & & `NPS_Dummy_nto_N_c1a' & & `NPS_Dummy_nto_N_c1b' & & `NPS_Dummy_nto_N_c2' \\
				tex \hline \hline 
				tex \end{tabular}
				tex {\footnotesize \onehalfspacing Notes: The table reports the doubly robust DID estimate of the effect of FBFs on each row's outcome variable for the cohort and subsample described in the column header. Columns 1-3 report the coefficients of interest for cohort-specific doubly robust DID regressions that exclude near-treated households (those between 4 and 5.5 hours from a FBF). The comparison is between households within the FBF service area (<4 hours travel time) and households far (>5.5 hours travel time) from a service area. Columns 4-6 report spillover coefficients for a doubly robust DID regression where the treated households are those located near but not within a FBF service area (between 4-5.5 hours) and the control households are those located far (> 5.5 hours away. For each comparison sample, we separately estimate a doubly robust DID estimate for each cohort as defined previously. The standard errors reported below each coefficient (in parentheses) are wild bootstrap standard errors clustered at the zone level. Stars reported alongside each coefficient represent the statistical significance: $* \left(p<0.1\right), ** \left(p<0.05\right), *** \left(p<0.01\right)$. The outcome's control group pre-period mean and standard deviation (in parentheses) are shown below the coefficient for each sample and cohort.}
				tex \end{threeparttable}
				tex \end{adjustbox}
				tex \end{table}
				texdoc close
					
					//Table S.20: Effect of FBFs on fertilizer adoption (top panel) and application rates (bottom panel) by cohorts as estimated using a continuous DID model.
						use FTF.dta, clear 
						gen T1_P1_cont= Treatment_1_continous*Post_1*(-1)
						gen T2_P2_cont= Treatment_2_continous*Post_2*(-1)
						gen T1_P2_cont= Treatment_1_continous*Post_2*(-1)

						loc outcome_all_datasets NPS_Dummy DAP_Dummy Fert_Dummy ///  // Fertilizer adoption - binary outcome   
							w_NPS_rate w_DAP_rate  w_App_rate

							*Cohort 1a 
							preserve 
							drop if year==3
							loc covariates_c1 land_area_sqm_10thou_y1 Female_y1 Litrate_y1
								foreach v in  `outcome_all_datasets' {
									wildbootstrap reg `v' Post_1 Treatment_1_continous T1_P1_cont `covariates_c1', cluster(A05) rseed(19810310) 
										estimates store `v'_c1a	
										}
							restore 

							*Cohort 1b
							preserve 
							drop if year==2
							loc covariates_c1 land_area_sqm_10thou_y1 Female_y1 Litrate_y1
								foreach v in  `outcome_all_datasets' {
									wildbootstrap reg `v' Post_2 Treatment_1_continous T1_P2_cont `covariates_c1', cluster(A05) rseed(19810310) 
										estimates store `v'_c1b	
										}
							restore 
							
							*Cohort two  
							preserve 
							drop if year==1
							drop if first_treat==2
							loc covariates_c2 land_area_sqm_10thou_y2 Female_y2 Litrate_y2
								foreach v in  `outcome_all_datasets' {
									wildbootstrap reg `v' Post_2 Treatment_2_continous T2_P2_cont `covariates_c2', cluster(A05) rseed(19810310) 
										estimates store `v'_c2
										}
							restore 	
							
						*Fertilizer adoption 
							*Cohort 1a and 1b
							esttab  NPS_Dummy_c1a DAP_Dummy_c1a Fert_Dummy_c1a ///
									NPS_Dummy_c1b DAP_Dummy_c1b Fert_Dummy_c1b ///
							using "Tables\Continuous_DID.rtf", ///
							nogap onecell  label se starl(* 0.10 ** 0.05 *** 0.01) mlabels(, dep) alignment(D{}{}{}) ///
							b(%09.3fc) order (T1_P1_cont T2_P2_cont) ///
							title(Table ##: Effect of fertilizer blending facilities on Fertilizer adoption: continous DID estimates: cohort 1a and 1b ) ///
							drop ( land_area_sqm_10thou* Female* Litrate* Treatment* Post*) replace   	
							
							*Cohort 2	
							esttab  NPS_Dummy_c2 DAP_Dummy_c2 Fert_Dummy_c2 ///
							using "Tables\Continuous_DID.rtf", ///
							nogap onecell  label se starl(* 0.10 ** 0.05 *** 0.01) mlabels(, dep) alignment(D{}{}{}) ///
							b(%09.3fc) order (T1_P1_cont T2_P2_cont) ///
							title(Table ##: Effect of fertilizer blending facilities on Fertilizer adoption: continous DID estimates: cohort 2 ) ///
							drop ( land_area_sqm_10thou* Female* Litrate* Treatment* Post*) append   
							

						*Application rates 
							*Cohort 1a and 1b
							esttab  w_NPS_rate_c1a w_DAP_rate_c1a  w_App_rate_c1a ///
									w_NPS_rate_c1b w_DAP_rate_c1b  w_App_rate_c1b ///
							using "Tables\Continuous_DID.rtf", ///
							nogap onecell  label se starl(* 0.10 ** 0.05 *** 0.01) mlabels(, dep) alignment(D{}{}{}) ///
							b(%09.3fc) order (T1_P1_cont T2_P2_cont) ///
							title(Table ##: Effect of fertilizer blending facilities on fertilizer application rate: continous DID estimates by cohort ) ///
							drop ( land_area_sqm_10thou* Female* Litrate* Treatment* Post*) append   
							
							*Cohort 2
							esttab  w_NPS_rate_c2 w_DAP_rate_c2  w_App_rate_c2 ///
							using "Tables\Continuous_DID.rtf", ///
							nogap onecell  label se starl(* 0.10 ** 0.05 *** 0.01) mlabels(, dep) alignment(D{}{}{}) ///
							b(%09.3fc) order (T1_P1_cont T2_P2_cont) ///
							title(Table ##: Effect of fertilizer blending facilities on fertilizer application rate: continous DID estimates by cohort ) ///
							drop ( land_area_sqm_10thou* Female* Litrate* Treatment* Post*) append 
				
				log close
				
				
